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ABSTRACT 


Under a contract with the National Aeronautics and Space Administration 
(NASA, Johnson Space Center in Houston) , Texas A§M University, through its 
Remote Sensing Center, is one of a small group of contractors, evaluating 
the feasibility of remotely sensing the status of soil moisture. 

At present, the most promising tedinology is the measurement of micro- 
wave emission by the earth's surface or the backs cat tering of microwave 
radiation emitted by an airborne source. In either case, the signals 
characterize a rather shallow layer, varioirsly estimated between 0.05 
and 0.25 m deep. Surface configuration, moisture content, or perhaps 
more accurately, the specific free energy of the water present, affect 
not only the signal strength, but also the depth of perception. 

To siqjport such efforts, a closely related stucfy of the moisture 
movement and the water balance of deeper layers is indispensable. This 
classical problem of agricultural and natural hydrology has always been 
made intractable by the lack of facts on the surface moisture regime. 

Thus, the two areas of research and application are currently merging. 
Evidence of this trend is the development of theories, in the form of 
numerical models, that are combinations of atmospheric and soil physics. 
These may proAade the foundation for models of microwave physics as well 
as a practical means to interpret radar signatures. 

As one step in that direction, we present here a conprehensive, yet 
fairly simple model of water disposition in a bare soil profile under the 
sequential inpact of rain storms and other atmospheric influences, as they 
occur from hour to hour. This model is intended mostly to support field 
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studies of soil moisture dynamics by our current team, to serve as a 
background for the microwave measurements, and, eventually, to serve as 
a point of departure for soil nwisture predictions or estimates based in 
part ipon airborne measurements. 

The main distinction of the current model is that it accounts not 
only for the moisture flow in the soil-atmosphere system, but also for 
the energy flow and, hence, calculates system tenperatures . Also, the 
model is of a dynamic nature, capable of supporting any required degree 
of resolution in time and space. 

This work is the. precursor of similar work for vegetated areas. It 
should be emphasized, however, that much critical testing of the sinple 
case is needed before the ccxrplexities of the hydrology of a vegetated 
surface can be related meaningfully to microwave observations. 

The present model is given in full detail, so as to invite its use by 
others and to make it possible to make adaptations, changes, and 
inprovements . The program, as listed, can be obtained from the authors 
in the form of punched cards or a cassette. 


IV 


INTRODUCTION 


A simulation model, is documented for calculating the water content and 
temperature profiles of a bare soil, from kiown initial conditions and a 
set of ordinary, time -dependent weather data, over a period of several days 
to several weeks. The model is dynamic, because the properties of the sys- 
tem are updated as the ten^rature and water content are changing in time. 

The present model was adapted from a set of algorithms devised for the 
study of dry mulching as a water conserving treatment in the dryland 
fanning areas of North Texas (Horton, 1977). In turn, the latter model 
was derived from a simulation of the concurrent flow of water and heat in 
soil proposed by Van Bavel and Hillel (197S, 1976) and of the infiltration 
and detention of rainfall suggested by Hillel et ail. (1975). 

The model provides a conprehensive method for the simultaneous solutim 
of the equation of continuity for water ahd heat in a soil system. The 
solution is obtained at frequent, fixed intervals and the moisture and 
tenperature profiles are printed when desired. The distinguishing 
characteristic of the model is that it does not assume a typical or average 
rate of evaporation, but that it, rather, generates the instantaneous rate, 
from the ambient weather and the momentary values of the soil moisture and 
temperatures. The evaporative flux is found by a combination method, that 
is, a combination of a surface energy balance and a model of the fluxes 
above and below that surface. 

The model is written in the Continuous System Modeling Program III 
(CSMP III) language (IBM, 1975), a specific numerical simulation language 
suitable for time- variant systems. The model was developed for execution 
on the AMDAHL 470V/6 computer operated by Texas A§M University at 



College Station, Texas. This system is operated with the IBM S-370 
operating system. 

A basic knowledge of physics and FORTRAN is desirable to understand 
the model. To assist the user, a glossary of variable names with their 
units has been provided as Appendix A, and a complete listing of the 
program is presented in Appendix B. Furthermore, since CSMP III is used 
throughout, a brief explanation will be given when C31P III statements 
are used. 
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IIJPUTS AND INITIAL VALUES FOR THE Ca-IPUTER SMULATION 

The inputs to calculate the heat and water balance of a soil system are 
divided in two parts: constant and variable inputs. 

Constant inputs to the model refer to the hydraulic cliaracteristics of 
the soil and the relationship between albedo (shortwave reflectance) and 
volumetric water content of the soil surface. The hydraulic cdiaracteristics 
of the soil are the functions relating pressure potential and hydraulic 
conductivity to volumetric water content. 

Variable inputs to the model refer to the time-dependeht weather variables. 
These are: 

1. daily global radiation, based on either daily total or hourly data, 

2. daily maximum and minimum air temperatures and their corresponding 
relative humidities, or, hourly tenperature ahd dewpoint data (2.0 m height) , 

3. average daily wind speed, or hourly data if available (2.0 m height) , 
and 

4. amount and duration of precipitation. 

Initial values that must be known are the initial soil moisture and 
heat content of the soil as a function of depth. These are obtained as 
follows : 

1. initial heat content are calculated from an initial measured or 
estimated soil temperature profile, and 

2. initial water contents as measured in the field. ' 

The way in which the constant and time -dependent inputs, and the initial 
values are manipulated in the model will be considered in detail in section 
IV of this document. 
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DESCRIPTION OF THE SOIL 


This simulation model is being applied to describe the heat and water 
balance of a soil classified in the Norwood series (Mixed (Calcareous) , 
Theimic Typic Udifluvents) . Specifically, the hydraulic characteristics in 
the model were obtained from three unpublished theses from Texas A§M Uni- 
versity (Saffaf, 1966; Marek, 1977; and Hunphreys, 1979). These reports 
deal with the field determination of the hydraulic conductivity of the soil 
and with the hydraulic characterisitics of irrigation furrows. The 
experimental work was done on the Agronomy farm of Texas A5M University, 
on a Norwood soil. 

in terms of its hydraulic characteristics, the soil profile being 
used in this simulation can be divided into two horizons (Hunphreys, 1979). 
The characteristics of each horizon are listed below: 

A. Horizon number 1: 


a. 

Depth 

0.0 - 0.20 m 

b. 

Texture 

silty clay loam 

c. . 

Average dry bulk density 

li.54 g/cn? 

d. 

Saturated hydraulic conductivity 

5.0 X 10 ^ m/s 

e. 

Relationship of pressure potential versus volumetric water 
content as plotted in Figure 1, and 

f. Calculated hydraulic conductivity versus volumetric water 
content as plotted in Figure 2. 

Horizon number 2: 

a. 

Depth 

0.20 - 1.20 m 

b. 

Texture 

silt loam 

c. 

Average dry bulk density 

1.46 g/cm^ 
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Figure 1 


Relationship between the Log-- of pressure potential 
(-m of water colunn) and voluifetric water content for 
the surface horizon, 0.0 - 0.20 m, of the Norwood soil 
(FUNCTION TVSPl) . 






LOGIO COND 


CO. 
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Relationship between the calculated hydraulic conductivity 
in m/s and the volumetric water content of the surface 
horizon, 0.0 - 0.20 m, of the Norwood soil (FUNCTION TVSCl) 






d. Average saturated hydraulic 6.0 x 10 ^ m/s 

conductivity 

e. Relationship of pressure potential versus volumetric water 
content as plotted in Figure 3, and 

f. Calculated hydraulic conductivity versus volumetric water 
content as plotted in Figure 4. 

The hydraulic conductivity as a function of volumetric water content 
for both horizons was calculated by the method of Jackson (Jackson, 1972). 
The WATFIV algorithm used for this calculation is included in Appendix C. 
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Figure 3. Relationship between the Log-» of pressure potential 
(-m of water column) and voltynetric water content for 
the sub-soil horizon, 0.20 - 1.20 m, of the Norwood 
soil (FUNCTION TVSP2) . 
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Figure 4. Relationship between calculated hydraulic conductivity 
in m/s and volumetric water craitent for the sub-soil 
horizon, 0.20 - 1.20 m, of the Norwood soil (FUNCTION TVSC2) 





DESCRIPTION OF THE NDDEL 


Introduction to CSMP III: 

Program statements in CSMP III can be divided into three categories: 
data, structure, and control statements. Data statements assign numeric 
values to the parameters, constants, and initial conditions of the system. 
Structure statements define the functional relationship between the 
variables of the model. Control statements deal with the duration of the 
simulation, the integration increment, type of output, and with the trans- 
lation and execution of the program. The user does not need to specify 
the category of each statement. 

The structure of the program can be divided into three sections: 

INITIAL, DYNAMIC, and TER^^INAL sections. The major characteristics of 
these translation control statements are the following: 

INITIAL Section 

1. Operations specified in this section are executed only at the 
onset of the simulation, 

2. it contains the equations that define the invariable geometry 
of the system, and 

3. it contains values and tables for specified parameters, and 
provides the initial values of specified variables. 

DYNAMIC Section 

1. It contains the equations that are needed to update tlie system 
at every time interval, and 

2. it uses an iterative procedure for the solutirai of inplicit 
functions of certain variables. 
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TERMINAL Section 


1. it specifies the finish time for simulation, 

2. it specifies the time interval for output and the variables 
to be printed in the output, and 

3. it specifies the method of integration that is to -be used 
and the integration time interval. 

Description 

The program will be described in the sequence given in the listing 
in Appendix B. To assist the user with the identification of variables 
refer to Appendix A where a glossary of the variables with their 
corresponding units are given. The International System of Units (SI) is 
used throughout the program with the exception of water potential values , 
which are given in terms of m water column. In principle, the water 
potential should be specified in'kPa, but the approximation that 
1 kPa = 0.1 m water column is sufficiently accurate and simplifies the 
dimensions of the units throughout. For the purpose of the descripticai, the 
program is divided into five parts: 

JOB CONTROL LANGUAGE (JCL) (lines 4-10) 

The JCL statements are the ones used at the Texas A^M University 
computer installation (AMDAHL 4 70V/ 6, IBM S370 system). The user should 
consult the instructions of the local installation and make the necessary 
changes. Note that the program requires at least 128K bytes of memory. 
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TITLE, NEMORY ORGANIZATION, AND ALLOCATION (lines 16-30) 

Lines 16-19 are TITLE statements. This is a CSMP III label statement 
and allows the user to specify a heading that will sppear at the top of 
the first page of printed output. Lines 14 and 15 are CSMP III comment 
statements, and can be identified by an asterisk (*) in column 1. Lines 
21-30 are translation control statements to organize the memory and to 
initialize arrays. Lines 21-23 are the CSMP III STORAGE statement, and 
they represent subscripted variables with the appropriate number of 
storage locations contained within the parentheses. Lines 24-27 are 
DIMENSION statements, and are handled as in FORTRAN with the exceptiem 
that a virgule (/) must appear in column 1. The virgule indicates that 
the DIMENSION instruction is a FORTRAN specification statement. Lines 27 
and 28 are EQUIVALENCE statements and are treated in the same manner as 
the DIMENSION statements. The EC^IVALENCE statement allows that the 
variables within the parentheses be assigned to the same storage locations 
in the memoiy, that is, the variables are synonymous. Line 30 are 
variables specified as integers with the statement FIXED. 

INITIAL Section (lines 36-210) 

The INITIAL section begins with the lines (the numbers to the left 
refer to the corresponding number of tJie listing in Appendix B) : 

36. INITIAL 

38. NOSORT 

The NOSORT is a CSMP III translation control statement and it means that 
the subsequent statements are to be executed in the order in which they 


12 


appear. The INITIAL section will be described in paragraphs as indicated 
by the comment cards in the listing (Appendix B). 


*** 1) DEFINITION OF PARAMETERS (lines 40-50) 

Lines 40 and 42-50 are data control statements. The PARAMETER statement 
(lines 40-50) is used to assign numerical values to the following variables 
(the number to the right in parenthesis refers to the corresponding number 
of the variable as given in the glossary in Appendix A): 


40. 

PARAMETER 

NL 

s 

13 

C68) 

42. 

PARAMETER 

KONDS 

ss 

4.20 

(61) 

43. 

PARAMETER 

KNOEW 

= 

0.57 

(63) 

44. 

PARAMETER 

KONDA 

= 

0.025 

(60) 

45. 

PARAMETER 

VHCAPS 

= 

1.925E-06 

(102) 

46. 

PARAMETER 

VHCAPW 

= 

4.186E06 

(103) 

47. 

PARAMETER 

SI04A 

= 

5.67E-08 

(82) 

48. 

PARAMETER 

ZO 

s 

0.01 

(108) 

49. 

PARAMETER 

SATCON 

s= 

0.60E-06 

(80) 

50. 

PARAMETER 

PORSTY 

= 

0.50 

(71) 


The values for the heat conductivity of soil (KONDS) , water (KONDW) , 


heat conductivity by air (KDNDA) and volumetric heat capacity of the soil 


(VHCAPS) and of water (VHCAPW) were obtained as suggested by De Vries (1966) . 
The value for the total porosity (PORSTY) of the soil was calculated 


from the average dry bulk density of the soil layers i.e, 

_ , n Average dry bulk density 
PORSTY - 1.0 pirticle dih^rty 


*** 2) DAILY COUNTERS (lines 52-54) 


Lines 52-54 initialize day counters. Which are used to keep track of 


the daily input data that is read with lines 56-61 in the third paragraph. 


*** 3) READ INPUT DATA (lines 56-61) 

The input data is stored in a two dimensional array (WINPUT) , and is 
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read in with a FORTRAN READ statement (line 57) . 


*** 4) SPECIFICATION OF THE GEOMETRY OF THE SYSTEM, INITIAI. WATER CONTENT 
AND TEMPERATURE (lines 64-68) 

In CSMP III, the TABLE statement is used to assign values to the 
subscripted variables listed on a STORAGE card. Thus, lines 64, 66, and 
67 are used to initialize the values for the thickness of the compartment 
(TCCM) , initial volumetric water content (ITHETA) , and initial soil 
tenperatures (ITEMP). Since the nun4)er of layers (NL) is 13, one value 
for each soil layer must be specified with the TABLE statement. Three 
periods (...) are used as continuation to following lines. 

I 

*** 5) CALCULATIONS OF DISTANCE AND DEPTH (lines 70-75) 

The depth (DEPTH(I)) and the distance (DIST(I)) of each soil layer is 
calculated in lines 70-75. The depth of each layer is the vertical 
distance between its midpoint and the soil surface, and the distance 
between midpoint of adjacent layers is given by DIST(I) (lines 72-75). 

The depth of the first layer (DEPTH(i)) is half its thickness (line 70) , 
and the distance of the first layer (DIST(l)) is set equal to its depth 
(line 7l) . 

70. DEPTH (1) = 0.5*TCOM(1) 

71. DIST(l) = DEPTH(l) 

72. DO 20 I = 2,NL 

73. DIST(I) = 0.5*(TCOM(I-1)+TCOM(I)) (21) 

74. DEPTH(I) = DEPTH(I-1)+0.5*(TCOM(I-1)+TCOM(I)) (17) 

75. 20 CONTINUE 

Figure 5 shows the geometry of the system and syn4)ols used in the 
program, and Table 1 lists the corresponding values of TOOM, DIST, and 
DEPTH. 
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Figure 5. Geometry of the system and symbols used in the program. 
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TABI^ 1. Invariable geometry of the soil system. 


Layer N 

(I) 

Thickness of 
conpartment 

TCOM(I) (m) 

Distance 
DIST(I) Cm) 

Depth 

DEPTH (I) Cm) 

1 

0.01 

0.005 

0.005 

2 

0.02 

0.015 

0.020 

3 

0.02 

0.020 

0.040 

4 

0.03 

0.025 

0.065 

5 

0.03 

0.030 

0.095 

6 

0.04 

0.03B 

0.130 

7 

0.05 

0.045 

0.175 

8 

0.05 

0.050 

0.225 

9 

0.10 

0.075 

0.300 

10 

0.10 

0.100 

0.400 

11 

0.15 

0.125 

0.525 

12 

0.30 

0.225 

0.750 

13 

0.30 

0.300 

1.050 
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*** 6) CALCULATIONS OF INITIAL WATER AND TEMPERATURE CONDITIONS (lines 79-87) 
From the parameters given, the following initial conditions for water 
and temperature in the soil system are calculated: 

1) the initial volumetric heat capacity of each soil layer (IVHCAP(I)) 
is calculated from the soil porosity and its water content as suggested by 
De Vries (1966) . 

81. IVHCAP(I) = VHCAPW*I’mETA(I) + (1.0-PORSTY)*VHCAPS (54) 

2) the initial water contmt (IWATER) for the entire soil profile is 
calculated by, 

82. IWATER = IWATER + T00M(I)*I'mETA(I) (57) 

3) the initial volume of heat of each soil layer, (IVOLH (I)) is 
calculated by, 

83. IVOLH(I) » ITEMP(I)*TCOM(I)*IVHCAP(I) (55) 

4) the net difference between influx and outflux of soil water (NFLUX) 
and inflow and outflow of heat (NFLOW) is set equal to zero respectively by, 

84. NFLUX(I) =0.0 (67) 

85. NFLOW(I) =0.0 (66) 

5) and the initial volume of water in each soil layer (IVDLW(I)) is 

calculated by , 

86. IVOLW(I) = nHETA(I)*T00M(I) (56) 

In CSMP III, a graphical relationship between pairs of x - y coordinates 

can be represented by the data control statement FUNCTION. Thus, lines 
93, 116, 139, 164, 188, and 198 represent the following relationships: 

1) FUNCTION TVSPl: Volumetric water content versus pressure potential. 

These data characterize the first soil horizon as given in Figure 1. 

2) FUNCTION TVSCl: Volumetric water content versus calculated 

hydraulic conductivity. These data characterize the first soil horizon as 
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given in Figure 2. 

3) FUNCTION TVSP2: Volumetric water content versus pressure potential. 

These data characterize the second soil horizon as given in figure 3. 

4) FUNCTION TVSC2: Volumetric water content versus calculated 

hydraulic conductivity. These data characterize the second soil horizon 
as given in figure 4. 

5) FUNCTION TEVSKO: Soil tenperature versus heat conductivity by 

water vapor. This relationship is plotted in Figure 6. 

6) FUNCTION TIVSAL: Volumetric water content of the first soil layer 
Versus albedo (shortwave reflectance). This relationship is plotted in 
Figure 7. 

*** 7) TABLE OF THE GEOMETRY OF THE SYSTEM AND ITS INITIAL STATE 
(lines 206-210) 

The INITIAL section ends with an output of a table of the geometry 
of the system as given in Table 1. 

DYNAMIC Section (lines 216-428) 

The equations which update the system at every time interval and to 
perform integrations are given in the DYNAMIC section. The Di’N/XMIC 
section will be described in paragraphs as indicated by the comment 
cards in the listing (Appendix B) . 

The DYNAMIC section begins with the lines , 

' 216. DYNAMIC 

218. NOSORT 
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Albedo versus volumetric water content of first soil layer 
CFUNCTION TIVSAL) . 
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As explained before, the NOSORT statement indicates that the statements 
are executed in the order in which they are listed. 

The paragraphs are: 

*** 1) DEFINITION OF TIME RELATED VARIABLES (lines 221-230) 

This paragraph defines the time related variables: time in hours 

(HTIME) and standard time of the day in hours (STIME). Also defined is 
an inpulse generator (IMPULS). 

223. Y = IMPULSC 86400., 86400.) 

which increments the day index (DNUM) , which is used to reference WINPUT. 

The first value in parentheses in tlie IMPULS statement is the time in 
seconds to the first pulse, and the second value is the time interval in 
seconds between pulses. 

The variables IMJM, DNUMl, and DNUM2 are daily counters that are used 
to keep track of the input data read frran the arr^ WINPUT where the data 
are stored. The method used to determine to vhat Julian day number (JDNUM) 
the input data corresponds is given by line 230. 

230. JDNUM = WINPUTi;i, DNUM) (58) 

*** 2) CALCULATION OF HYDRAULIC CHARACTERISTICS OF THE FIRST HORIZON 
(lines 232-237) 

The volumetric water content Cn-ETA(I)), hydraulic conductivity (COND(I)), 
pressure potential (PPOT(I)), and hydraulic potential (HPOT(I)), for the 
first 8 layers of the soil profile (horizon #1) are calculated from the 
following equations: 

1) 'IHETA(I) of each layer is obtained from the ratio of the 
volume of water (VOLW(I)) to the layer volume per unit area: 
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233. 


THETA(I) = VOLW(I)/TOOM(I) 


(93) 


2) CDND(I) at the center of each layer is obtained by linear 
interpolation of the values in the FUNCTION TVSCl using the water content 
(THETACI)) of each layer. 

234. COND(I) = AFGEN (TVSCl, TOETA(I)) (7) 

AFGEN (Arbitrary Function ffiNerator ) is a CSMP III statement that allows 
linear interpolation from a relationship specified with a FUNCTION 
statement. 

3) PPOT(I) is obtained by linearly interpolating values in the 
FUNCTION TVSPl using THETA(I) of each layer. 

235. PPOT(I) = AFGEN (TVSPl, THETAfl)). (72) 

4) HPOT(I) is the sum of PPOT(I) and the elevation (-EEPTH(I)) 
of each layer. 

236. HPOT(I) = PPOT(I) - DEPTH(I) (45) 

*** 3) CALCULATION OF HYDRAULIC CHARACTERISTICS OF SECOND HORIZON 

(lines 239-244) 

The same calculations outlined in paragraph #2 are performed for the 
bottom five layers of the soil profile, that is, the second horizcn. In 
this case, FUNCTIONS TVSP2 and TVSC2 are xised with the AFGEN statements . 

*** 4) CALCULATION OF THERMAL PROPERTIES (lines 246-253) 

In this paragraph, the Volumetric heat capacity (VHCAP(I)), temperature 
(TEMP(I)), and thermal conductivity (KOND(I)) are calculated for the 13 
layers in the soil profile. Tlie equations used are: 

1) VHCAl’(I) is calculated from the porosity (POllSTY) and the water 
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content (THETA(I)) as suggested by De Vries (1966) . Soil air is not 
considered in the calculation. 

247. VHCAP(I) = VHCAPW*THETA(I) + (1.0 - PORSTY) *VHCAPS (101) 

2) TEMP(I) of each soil layer is calculated by dividing the 
volumetric heat content per unit area (VOLH(I)) by the product of 

volumetric heat capacity (VHVAP(I)) and layer thickness (T0C»4(I)). 

248. TEMP(I) = VOLH(I)/(VHCAP(I)*TCOM(I)) (91) 

3) The contribution of the water vapor phase to the thermal 
conductivity (KONDV) is obtained by linear interpolation of values in 
FUNCTION TCVSKO and TEMP (I) for each layer. 

249. KONDV =AFGEN(TEVSKO, TEMP (I)) (62) 

4) KOND(I) of each soil layer is found by the formula suggested by 
De Vries (1966) using the values of KONDS, KCNDA, and KONDW assigned in 
the INITIAL section. 

250. KOND(I) = ((1. -PORSTY) *KONDS*0. 4 + THETA(I)*KONDW +... (59) 

(PORSTY - THETA(I))*1.4*(KONDA + KNODV))/... 

( (1. -PORSTY) *. 4+THETA(I) + (PORSTY-'mETA(I) ) *1 . 4) 

*** 5) CALCULATIONS OF AVERAGE CONDUCTIVITIES (lines 255-260) 

The average thermal conductivity (AVKOND) and average hydraulic 
conductivity (AVCOND) for .transport between adjacent layers is taken to 
be the average of layer conductivities, weighted according to their relative 
thickness. 

For the flow of heat, AVKOND(i) is calculated from: 

256. AVKOND(I) = (TO)M(I-1)+TCOM(I))/(TCOM(I-1)/KOND(I-1)+TCON1(I)/ (4) 

KOND(D) 

For the flux of liquid water, AVCOND(I) is calculated from: 
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258. AVC0ND(I) = (C0ND(I-1)’'TG0M(I-1)+C0NDCI)*TCXDM(I))/(T00M(I-1) (3) 

+TOOM(I)) 

*** 6) SPECIFICATIONS OF BOUNDARY CONDITIONS FOR FLOW OF HEAT AND 
• FUJX OF WATER (lines 263-264) 

Tlie boundary conditions at the bottan layer (NLL) of the soil profile 
are defined by the following: 

1) The flux of water at the bottom boundary is taken to be equal to 
the hydraulic conductivity of the bottom layer. That is, the flux of water 
is driven by mit hydraulic potential gradient 

263. FLUX (NLL) = COND(NL) (39) 

2) The flow of heat at the bottom boundary is calculated by Fourier’s 
law assuming that the temperature at 1.20 m depth remains constant, and 
that this ten^ierature is set equal to ITEMP(NL). 

264. FLOW(NL) = (TEMP(NL) - ITEMP(NL) * KOND(NL/TCOM(NL)/2.0) (37) 

*** 7) CALCULATION OF FLOW OF HEAT AND FLUX OF WATER (lines 266-269) 

The flow of heat between layers is calculated by Fourier’s law. 

267. FLOW(I) = (TEMP(I-l) - .TEMP(I)) * AVKOND(I)/DIST(I) (36) 

The flux of water between layers is calculated by Darcy's law. 

268. FLUX(I) = (HPOT(I-l)-HPOT(I)) * AVCOND(I)/DIST(I) (33) 

Note that flow of heat and flux of water is calculated for all layers 

except the top one . 

*** 8) USE DAILY RAINFALL DATA, Af4D CALCULATE DAILY RAINFALL DISTRIBUTION 
(lines 272-285) 

The daily rainfall distribution is assumed to follow the 
pattern illustrated in Figure 8. The necessary inputs to calculate the 
rainfall distribution are given by. 
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RflIN (MM/HR) 



TIME OF DRY (HOURS) 


Figure 8. Pattern of daily rainfall distribution. For definition 
of symbols see text. 
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272. 

BEGIN = WINPUT(9, DNUM) 

(6) 

273. 

END = WINPUT(10, DNUM) 

(32) 

274. 

RFT = WINPUT(11, DNUM) 

(75) 


The logical branching used to calculate the rainfall distribution as 
a function of time is illustrated by the flow chart of Figure 9. In the 
program this is calculated by lines 272-285. 


*** 9) USE DAY LENGTH AND CALQJLATE DAILY DISTRIBUTION OF GLOBAL 
RADIATION Clines 288-292) 

The daily distribution of global radiation (GR) is spread out over the 
daylength tUL) period using a sine function. An exaiiple of this distribution 
is given in Figure 10. 

290. GR = 436.33*WINPirr(3, DNUM)/DL * SINC(STIME -12. + DL/2.) (41) 

*3.141/DL 

Note: 436.33 is a sinplification of (24x10^86400) x(it/2) 

***10) CALCULATION OF ALBEDO (lines 294-295) 

Albedo (ALB) is obtained by linear interpolation of values in TABLE 
T1V3AL and the volumetric water content of the first soil layer (Tl) . 

294. Tl = THETA(l) 

295. ALB = AFGEN (TIVSAL, Tl) 

***11) CALCULATION OF WINDSPEED AND BOUNDARY LAYER RESISTANCE (lines 297-305) 
Values of the average windspeed (SA) are set at noon (line 297) and 
linear interpolation produces values at other times (lines 299-302). The 
boundary layer resistance (RA) is calculated as the quotient of the 
logarithm of 2.0 divided by the roughness factor (ZO) squared and the 


(85) 

( 2 ) 


26 




Figure 9. Flow chart showing the logical branching used to calculate 
the pattern of rainfall distribution. 
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^.00 6.00 12.00 18.00 24.00 

TIME, .HRS 


Figure 10. Example of daily distribution of global 2 radiation (GR) 
over tine. (Total irradiance 10.0 MJ/m , day length = 
10.7 hours). 




product of windspeed at 2.0 meters and a stability factor (Sellers, 1965). 

305. RA » (ALOG(2.0/Z0)**2.0)/(0.16*SA) 

( 73 ) 

***12) CALCULATION OF DEWPOINT TEMPERATURE AND ABSOLITTE HUMIDITY 
(lines 307-314) 

Dewpoints at time of minimum air tenperature (Dl’MIN) are set at 5; 00 hours 
standard time and dewpoints at time of maximum air tenperature (DfWAX) are set 
at 15:00 hours standard time allowing linear interpolation of dewpoints 
(DPTC) at other times (line 311 and 313). 

The absolute humidity of the air (HA) is calculated from the 
equation (Murray, 1967) , and is plotted in Figure 11. 

314. HA = 1.323*EXP(17.27*DPTC/(237.3*DPTC))/(273.16+TAC) (42) 

***13) CALCULATION OF TEMPERATURE OF TIE AIR AND VOLUMETRIC HEAT 
CAPACITY OF HE AIR (lines 316-324) 

Minimum air temperature values (TAMIN) are set at 5:00 hours standard 
time and maximum air tenperatures(TAMAX) are set at 15:00 hours standard 
time allaving linear interpolation of air temperatures (TAG) at other 
times. 

Tlie volumetric heat capacity of the air (SH) , as a function of air 
temperature in degrees Kelvin (TAK) is calculated frort the following 
equation, and is plotted in Figure 12. 

324. SH » (1154..8*303.16)/(TAK) (81) 

***14) CALCULATION OF SKY IRRADIANCE (line 326) 

The sky irradiance (SKL) is found from a formula suggested by 
Brunt (1932), as modified by Sellers (1965). 
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TEMPERRTURE, C 


Figure 11. Relationship between absolute humidity of the air (HA) 
and tenperature as calculated with Murray's equation. 









326. SKL = (SI(MA*TAK**4)*(0.605 + 0. 039*SQRT(410. *HA) ) (83) 

***15) IMPLICIT CALCULATION OF THE SOIL SURFACE TEMPERATURE (lines 328-335) 
The tenperature of the soil surface (TSC) is calculated by an iirplicit 
function (line 328) using the air tenperature (TAC) as a first estimate and 


calculating the energy balance at the surface: 

1) Sensible heat flux (A) 

329. A = (TSC-TAC)*SH/RA ( 1 ) 

2) Murray's equation (Murr^, 1967) to calculate the saturated 
humidity (HO) at the soil surface. 

330. HO = 1.323*EXP(17.27*TSC/(273.3+TSC))/(273.16.TSC) (44) 

I 

3) Absolute humidity is calculated as sjiggested by Van Bavel and 
Hillel (1976). 

33. HO = HO*EXP(PPOT(l)/(46.97*(TSC+273.16)) (44) 

4) The evaporation (EV) is the quotient of the difference in humidities 
of surface and atmosphere and the boundary layer resistance. 

332. EV = (HO -HA)/(RA*1000.) (33) 

5) The conduction into the soil (S) is calculated as the difference 
in net radiation and the sum of sensible and latent heat fluxes. 

333. S = GR*(1.-ALB) + SKL-SIGMA(TSC+273. 16) ** 4-A-LH*EV (78) 

6) The final tenperature of the soil surface (FTSC) as calculated in 
the implicit loop is given by Fourier's law for conduction. 

335. FTSC = TEMP(l)+S*DEPTH(l)/KOND(l) (40) 


***16) CALCULATIONS OF EVAPORATION AND NET RADIATION (lines 337-344) 

The flow of heat into the center of layer one (FLOW(l)) is calculated 
as the product of heat conductivity (KOND(l)) and the differences in 
surface (TSC) and layer tenperatures (TEMP(l)) divided by 
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distance (DIST(l)). 

337. FLCHV(l) = (TSC-TEMP(1))*K0ND(1)/DIST(1). 

The saturated humidity of the air at the soil surface (HS) is 
calculated with htirray's equation. 

338. HS = 1.323*EXP(17.27*TSC/(237.3+TSC))/(273.16+TSC) (46) 

ITie absolute humidity (HS) is calculated by^ 

340. HS = HS*EXP(PP0T(1))/(46.97*(TSC+273.16)) (46) 

Evaporation (EVAP) is the quotient of the difference in humidities of 
soil surface and the atmosphere and the boundary layer resistance. 

342. EVAP = (HS-HA)/(RA*1000.0) (34) 

I 

The latent heat of vaporization (LH) as a function of the soil surface 
tenperature is given by (Forsythe , 1964) . 

343. LH - 2.94963E09 - 2.247E06*TSC (64) 

and in as a function of tenperature in plotted in Figure 13. 

Net radiation (NR) at the soil surface is calculated by difference 
from the energy balance equation. 

344. NR = FLOW(l) + (TSC-TAC)*SH/RA+LH*EVAP (70) 

A diagram illustrating the energy and water flux in the top two layers 

of the soil system is given in Figure 14. 

***17) CALCULATION OF DETAIN, INFILTRATION, AND INCAP (lines 346-359) 

The amount of water that remains on the soil surface (DETAIN) is 
defined as the integral between the initial value and the difference 
between the rainfall rate (RAIN) and tlie infiltration rate (INFILT) . In 
CSMP ill, this is acconplished by an INTGRL statement. 

346. DETAIN = INTGRL (0.0, RAIN - INFILT) (18) 

The infiltration capacity (INCAP) is the Darcian flux to the center 
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Figure 14. Schematic diagram of the energy and mass 
two layers of the soil system. 
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of the first layer from the saturated soil surface, at which the pressure 
potential is assigned a value of zero. 

347. INCAP = (0.0 - HPOT(1))*0.5*(SATCON+COND(1))/DIST(1) (49) 

Since the model assumes no runoff, only two possibilities exist with 
rainfall, and these are: 

1) IVhen the rainfall rate does not exceed the infiltration capacity 
and there is no detention on the soil surface, then the rainfall rate 
controls the infiltration rate. 

358. IF(RAIN.LT.INCAP.AND.DETAIN.LE.O) INFILT=RAIN 

2) When the rainfall rate exceeds the infiltration capacity of the 
soil ahd water is detained on the soil surface, the infiltration capacity 
controls the infiltration rate. 

357. INFILT = INCAP (51) 

The logical branching of the above calculations is illustrated by the 
flow chart in Figure 15. 

***18) CALCULATION OF NET FLOW OF HEAT AND NET FLUX OF WATER (lines 361-365) 
The flux of water into the middle of the soil surface layer is equal 
to the difference in the rate of infiltration and the rate of 
evaporation . 

361. FLUX(l) = INFILT- EVAP 

Tlie heat flow into the bare soil surface (FLOW(l)) was previoasly 
defined as, 

337. FLOW(l) = (TSC-TEMP(l))*KOND(l)/DIST(l) 

The flux of water (FLUX(I)) and flow of heat FIXHV(I)) for the rest of 
the layers must obey the continuity equation that is, 

Nnf)W(I) = FLOW(I) - FLOW(I-l) (66) 
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Figure 15. Flowchart showing the logical branching used to calculate 
INFILT and DETAIN. 


3 



NFLUX(I) = FLUX(I) - FLUX(I+1) 


(67) 


***19) INTEGRATION OF VOI^ETRIC HEAT CONTENT AND VOLUNETRIC WATER 
OONIENT (lines 368-369) 

The 13 integrations to keep track of the volumetric heat contents 
(VOLH(I)) and volumetric water oxitents (VOLW(l)) arc carried out by the 
CSMP III INTGRL function. 

368. VOLHl = INTGRL (IVOLHl,NFLOitfl. 13) (105) 

369. VOLWl »= INTGRL(IVOLWl,NFLUX1.13) (106) 

The third argument of the integral function indicates that there are 

13 integrals to keep track of the volumetric heat content and water 
content of the 13 layers. The heat contents are stored in an array VOLH(I) 
and the water contents are stored in an array VOLW(I). Note that these 
arrays are named in the INITIAL section (lines 28 and 29). The first 
argument of the integral function states that the initial value of the 
volumetric heat content and volumetric water content is given by an array 
IVOUI(I) and IVOLW(I), respectively. The second argument states that the flow 
and flux rate into the integral is given by the array NFLOW(I) for heat and 
NFLUX(I) for water. 

***20) CALCULATION OF CU^ULATIVE RAIN, INFILTRATION, EVAPORATION, AND 
Dl^AlNAGE (lines 372-375) 

Cumulative rain (CUMRN) , cumulative infiltration (CIMINF) , cumulative 
evaporation (CUNEVP) , and cumulative drainage (CUMURN) are calculated 
with the eSMP III INTGRL function.. 

372. CUMRN = INTGRL(0.0, RAIN) (11) 

373. CUMINF = INTGRL (0.0, INFILT) (10) 


38 



CUMEVP = INTGRUO.O, EVAP) (9) 

375. CUMDRN = INTGRLCO.O, FLXNLL) (8) 

***21) CALCULATION OF TOTAL WATER FOR DIFFERENT LAYERS (lines 377-392) 

'.The integrated volume of water in the profile , fron a given depth to 
the bottom of . the profile, is the sum of the volume of water in each layer 


starting at the specified depth. • 

1) From 0.0 - 1.20 m (CUNWTR) (1=1, NL) 

382. CUNWTR = CUNWTR+VOLW(I) (12) 

2) From 0.13 - 1.20 m (CUNWTl) (I = 6,NL) 

385. CUMWTl = ClMVTl + VOLW(I) , (13) 

3) From 0.30 - 1.20 m (aJM^T2) (1=9 ,NL) 

388. CUNWT2 = OMYr2+VOLW(I) (14) 

4) From 0.75 - 1.20 m (CUNWT3) (1=12, NL) 

391. CUMWT3 = CUMWT3 + FLOW(I) 


***22) CALCULATION OF DAILY TOTALS (lines 394-404) 

The daily totals of infiltration (DINF) , rain (DRN) , evaporation 
(DEVP) , and drainage (DDRN) are calculated by difference from the 
accumulated totals and the values at the previous midnight. 


400. 

DINF = CUMINF 

- INF(DNUM - 

2) 

(20) 

401. 

DRN = CUMRN 

- RN(DNUM - 

2) 

(30) 

402. 

DEVP = CUMEVP 

- EVP(DNUM - 

2) 

(19) 

403. 

DDRN = CUMDRN 

- DRAIN(DNUM 

- 2) 

(UO 


Tlie daily totals are calculated rally at the end of a day and tlie 
calculations are controlled by the IMI’ULS statement. 

394. ZBHJS = IMPULS(86400.,86400.), 
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***23) CHECK WATER BALANCE (line 406) 

The water balance for the soil system (BALANS) is calculated as the 
difference in the sum of cumulative water in the profile (CUNTTWR) , 
cumulative ev^oration (CUMEVP) , and cumulative drainage (CUMDRN) and the 
sum of the initial water content (IWATER) and the cumulative infiltration 
(CUMINF) . 

406. BALANS = CIMVTR - IWATCR - CUMINF + CUMEVP + CUMDRN (5) 

The deviatim of BALANS from zero provides an indication of the error 
produced in the solution of the equations in the model. 

***24) OUTPUT OF DESIRED VARIABLE (lines 408-428) 

The last paragraph of the DYNAMIC section deals with output of 
calculated variables. Specifically, the following variables are printed 
at 8:00 and 16:00 hours, standard time: DEPTH(I), IHETA(I) , PPOT(I), 

FLUX(I), NFLUX(I), andTEMP(I) 

The output at 8:00 hours is controlled by the IMPULS statement. 

408. Z = IMPULS (28800. 0. , 86400.) 

and the output at 16.0 hours is controlled by' the IMPULS state- 
ment. 

42S. ZZ = IMPULS (57600. 0. , 86400.0) 

TERMIXAL Section (lines 455-466) 

The terminal section is mainly conposed of execution control state- 
ments. The statements that appear in the TERMINAL section are: 

‘ 1) TIMER FINTIM = 259200., PRDEL = 84600.0, DELT = 100.0 
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The TIMER card is used to specify the variables that control the run-time, 
print increment, and integration interval (step-size): 

a) FINTIM determines the value of TIME (independent variable) at 
which the run is terminated. FINTIM is set equal to the desired 
simulation time. 

b) PRDEL. This TIMER variable controls the increment for the output 
of the PRINT statement. 

c) BELT. This TIMER variable specifies the integration interval. 

2) PRINT. The PRINT card is used to specify the variables that will 
be printed at each specified interval (PRDEL) during the simulation. 

3) METHOD TRAPZ. The integration technique is specified by the use 
of the METHOD card with the appropriate CSMP III integration name. In 
this case, the selected method is the trapezoidal (TRAPZ) method that 
uses a fixed integration interval (BELT). 

4) END. The END card specifies the end of the program. 

The last segment of the program is a list of the input data stored in the 
two dimensional array WINPUT. It is specified between the lines INPUT 
and ENDINPUT. This array must match the FORMAT statement of line 58. 

The program ends with the lines 

464. STOP 

465. ENDJOB 

466. /*END 

The ENDJOB card must begin in the first column. 
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DESCRIPTION OF THE OUTPUT OF THE MODEL 

To provide the user with an example of the output of this model, 
the program described in section IV of this document was executed for a 
period of three d^s. 

The output of the model consists of two parts: 

1) The first part is generated in the DYNAMIC section (paragr^h 24) , 
and for each simulated day an output at 8:00 am and 4:00 pm is printed. 

The variables printed for each soil layer are: DEPTH, THETA, PPOT, 

FLUX, NFLUX, and TEMP. 

2) The second part of the output is generated with the PRINT card 

the TERMINAL section. The variables printed for midnight at the end of 
each day of the simulation period are: JDNUM, XDNUM, DRN, DINF, lEVP, 

DDRN, CUMVTR, CUNWTl, OMVT2, OMVT3, BALANS, CUMRN, CIJMDRN, CUMINF, 
CIJMEVP, and FLOW(14). 

It should be noted that for a fixed integration method, such as 
TRAPZ, the integration interval (DELT) necessary for execution of the 
model may have to be adapted to the input data. Our experience with 
this model, using TRAPZ, suggests a DELT of 100 seconds for input data 
with rain, and 200 seconds for input data with no rain. However, the 
user has a choice of other integration methods and should consult the 
CSMP III manual for proper use. 
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APPENDIX A - GLOSSARY OF VARIABLE NAMES 
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riTL3: GLOSSARY FOR HEAT AND HATER 3ALANCE MODEL 


NO TERM DEFINITION. UNITS 

1. A SENSTDLE HEAT FLUX INTO THE AIR 

2. ALB ALBEDO 

3. AVCOND (I) .. AVERAGE HYDRAULtC CONDUCTIVITY 

.. BETWEEN ADJACENT LAYERS .M/S 

U. AVKOND(I) .. AVERAGE THERMAL CONDUCTIVITY 

..BETWEEN ADJACENT LAYERS W/(M*C) 


5. BALANS WATER BALANCE M 

5. BEGIN BEGINNING OF RAINFALL PERIOD HOURS 


7. 

8 . 

9 . 

10 . 
11 . 
12 . 
13. 
1 «. 
15. 


COND (I) . . . . 

CUMDRN 

COMEVP 

CUMINF 

CUMRN 

CUMWTR 

CUMWT1 

C0KHT2 

CUMBT3 


HYDRAULIC CONDUCTIVITY OF LAYER (I). M/S 


CUMULATIVE DRAINAGE M 

CUMULATIVE EVAPORATION M 

CUMULATIVE INFILTRATION M 

CUMULATIVE RAIN M 

TOTAL WATER IN SOIL PROFILE M 

TOTAL WATER (1=6, NL) M 

TOTAL WATER (1 = 9, NL) M 

TOTAL WATER (1 = 12, NL) M 


16. 

17. 

18. 

19. 

20 . 
21 . 

22 . 

23. 

24. 

25. 

26. 

27. 

28. 

29. 

30. 
3 1 . 


DDRN DAILY DRAINAGE 

DEPTH (I) • • . V MKIMC A L DISTANCE BEIKEEK iilDPCINT.. 

...OF LAYER (I) AND THE SURFACE 

DETAIN DEPTH OF HATER ON THE SOIL SURFACE.. 

DEVAP DAILY EVAPORATION 

DINF DAILY INFILTRATION 

DIST (I) .... DISTANCE BETWEEN MIDPOINTS OF 


DL DAY LENGTH 

DNUM DAY COUNTER 

DNUMI DAY COUNTER 

DN0N2 DAY COUNTER 

DPMAX DEWPOINT TEMP. AT TIME OF TMAX 

DPKIN DEWPOINT TEMP. AT TIME OF TMIN 

DPTC DEWPOINT TEMPERATURE 

DRAIN DRAINAGE 

DRN DAILY PAIN 

DKSLOP. . . . . P AT HFALI, SLOPE BETWEEN MIDPOINT AND. 

END. 


N 


M 

M 

M 

H 

M 

HOURS 



MK/HR<‘*2 


32. 

33. 

34. 

35. 


END END OF RAINFALL 

EV ....EVAPORATION (DU 

EV AP EVAPORATION 

EVP EVAPORATION (DU 


PERIOD HOURS 

MMY VARIABLE) M/S 

M/S 

MHY VA RIABLE) . M/S 


36. FLOW (I) FLOW OF HEAT INTO LAYER (I) W/M**2 

37. FLOW (NLL) .. FLOW OF HEAT ACROSS THE LOWER 

, . BOUNDARY W/M**2 

38. FLUX (I) . . . . FLUX OF WATER INTO LAYER (I) M/S 


39. FLUX (NLL) .. FLUX OF HATER ACROSS THE LOWER 
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no. 


<4 1 . 


< 42 . 

43. 

44. 

45. 

46. 
4 7. 
48. 


49. 

50. 

51. 

52. 

53. 

54. 

55. 

56. 

57. 


58. 


59. 

60. 
61. 
62. 
63. 


64. 


65. 


66 . 

67. 

68 . 

69. 

70. 


71. 

72. 


73. 


FTSC 


BOtJHDftflY M/S 

FINAL TEMPERATURE OF SOIL SURFACE... 

AS CALCULATED IN THE IMPLICIT LOOP. .C 


GR GLOBAL RADIATION 


H/M**2 


HA 

.ABSOLUTE HUMIDITY OF THE AIR 

. . KG/M**3 

HEIGHT. . . . 

.RAINFALL HATE AT MIDPOINT 

. .HM/HOUR 

HD 

.SATURATED HUMIDITY AT THE SOIL.... 
.SURFACE AMO ABSOLUTE HUMIDITY 

• • 

. . KG/K**3 

HPOr (1) . .. 

•HYDRAULIC POTENTIAL OF SOIL LAYER. 

..M OF HATER 

HS 

- ABSOLUTE. HUKTDTTY 

. . KG/M**3 

HSO. ...... 

.ABSOLUTE HUMIDITY, DUMMY VARIABLE. 


HUMS 

.TIME 

. .HOURS 


INCAP 

INF 

INFILT 

ITEMP (I) . . . 
IFHETA (I) . . 

• • 

IVHCAP (I) . . 

IVOLH (I) .. . 
IVOLW (I) ., . 

« • • 

IHATER 


INFILTRATION CAPACITY M/S 

INFILTPATION M 

INFILTRATION RATE H/S 

INITIAL ""EMPERATURE OF LAYER (I)....C 
INITIAL VOLUMETRIC HATER CONTENT.... 

OF LAYER (I) ..M**3/H**3 

INITIAL VOLUMETRIC HEAT CAPACITY.... 

OF LAYER (I) J/(M**3*C) 

INITIAL AMOUNT OF HEAT IN LAYER (1).J/M**2 

INITIAL VOLUME OF WATER IN 

LAYER (I) M 

INITIAL TOTAL WATER CONTENT OF 

THE SOIL PROFILE M 


JDNUH JULIAN DAY NUMBER 


KOND(I) 

RDNDA 

KONDS 

KONDV 

KONDW 


THERMAL CONDUCTIVITY OF LAYER (I) . . . H/ (MfC) 

THERMAL CONDUCTIVITY OF AIR W/(M'»C) 

THERMAL CONDUCTIVITY OF SOIL. ....... U/ (M*C) 

THERMAL CONDUCTIVITY BY HATER VAPOR. K/ (H^C) . 

THERMAL CONDUCTIVITY BY HATER H/(M*C) 


LH LATENT HEAT OF VAP. OF HATER J/M«‘*3 


MDPNT .MIDPOINT OF RAINFALL PERIOD. HOURS 


NFLOH (I) . . .NET FLOW OF HEAT INTO LAYER (I) V/K**2 

NFLUX (I) .. .NET FLUX OF WATER INTO LAYER (I) M/S 


NL HUMBER OF LAYERS - 

NLL. ....... NL *■ 1 - 

NR NET RADIATION AT THE SOIL SURFACE. .. W/K ♦♦2 


PDRSTY. .... POROSITY OF THE SOIL - 

PPOT (I) .... PRESSURE POTENTIAL OF LAYER (I) M OF HATER 


RA BOUNDARY LAYER RESISTANCE S/M 
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7U. 

75. 

75. 

78. 


78. 

79. 

80. 

81. 

82. 

83. 

84. 


85. 

86 . 

87. 

88 . 

89. 

90. 

91. 

92. 

93. 

94. 

95. 

96. 

97. 

98. 

99. 


100 . 


101 . 
102 . 
10 3. 
104. 
1 05. 


106. 


107. 


RAIN RAINPAI.L RATP . . N/S 

RFr TOTAL RAINFALL DETWEEN BEGIN AND 

END MM 

RHS RELATIVE HOMIDITI OF THE SOIL SURFACE - 

RN RAIN (DUMMY VARIABLE) .....M 


S 

....CONDUCTION OF ENERGY INTO THE SOIL.. 


SURFACE 

. . . . W/M*»2 

SA 

. . . . WINDSPEED 

. . . . M/S 

SATCON. 

....SATURATED HYDRAULIC CONDUCTIVITY 

OF. 


....THE SUOFACr. HOPTV.ON 


SH 

....VOLUMETRIC HEAT CAPACITY OF THE 

AIR. J/ (M*#3*C) 

. . . .W/(H**2*K<‘*4) 

A « t 

. * . - ^T^PAN-ROT.T?. MAM rnN<\T AMT 

SKL. . . . 

. . . .SKY RADIANCE 

....W/M'**2 

STIME. . 


. . . . HOURS 


n WATER CONTENT OF FIRST LAYER 

(DUMMY VARIABLE) M**3/M**3 

lAC TEMPERATURE OF THE AIR C 

TAK TEMPERATURE OF TUB AIR K 

TAHAX MA)!TMUn AlP TEMPERATURE C 

TAMIH MINIMUM AJR TEMPERATURE .C 

TCOM (1) ... .THICKNESS OF LAYER (I) M 

TEMP (I) ... .TEMPERATURE OF LAYER(I) C 

TEV3KO SOIL TEMPERATURE VS. THERMAL 

CONDUCTIVITY BY WATER VAPOR C VS. W/(M’CC) 

THETA (I) ... VOLUMETRIC WATER CONTENT OF LAYER... 

... (I) • M**3/M**3 

T1VSAL VOLUMETRIC WATER CONTENT OF FIRST... 

SOIL LAYER VS ALBEDO 

TSC TEMPERATURE OF THE SOIL SURFACE C 

TVSC1 VOLUMETRIC WATER CONTENT VS. HYDRAULIC 

......CONDUCTIVITY FOR FIRST HORIZON 

TVSC2. ..... VOLUMETRIC WATER CONTENT VS. HYDRAULIC 

CONDUCTIVITY FOR SECOND HORIZON 

TVSP1 VOLUMETRIC WATER CONTENT VS. PRESSURE 

POTENTIAL FOP FIRST HORIZON 

TVSP2 VOLUMETRIC WATER CONTENT VS. PRESSURE 

POTENTIAL FOR SECOND HORIZON 


OPSLOP RAINFALL SLOPE BETWEEN BEGIN AND 

MIDPOINT 


MM/HR**2 


VHCAP (I) ... VOLUMETRIC HEAT CAPACITY OF LAYER ( I) . J/ ( M ♦ ♦ 3* C) 

VHCAPS VOLUMETRIC HEAT CAPACITY OF THE SO I L . . J/ ( M ♦ * 3 *C) 

VHCAPW VOLUMETRIC HEAT CAPACITY OF W ATER . . • . . 0/ ( K ♦ ♦ 3 ^C) 

VOLH (I) . . . . VOLUMBTR IC HEAT CONTENT OF LAYER (I)..J/H^*2 

VOLW(I) VOLUME OF WATER P3R UNIT AREA OF LAYER 

(I) M**3/M**2 


WINPUT ARRAY FOP WEATHER INPUT DATA 


ZD SURFACE ROUGHNESS COEFFICIENT 
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APPENDIX B - PROGRAM LISTING 



1. «««* A JOB CONTRO. LANGUAGE 

2 . ***♦ 

3. 

A, //SCSDRYJA JOB ( RO A2. A03A ,003 t 020. RL ) . ' VAN BAVEL - LASCANO* 

5. /ALEVEL 0 

6. /AJOBRARM R=128 

7. // EXEC CSMD3CLG 

e. //COMPRINT 00 DUMMY 

9. //SYSPRINT 00 DUMMY 

JO. //SYSIN DO * ' 


11. **** 

12 . **** 

13. **** TITLE. MEMORY ORGANIZATION AND ALLOCATION 

14. *44* 

15. 4*4* 

16. TITLE WATER AND HEAT BALANCE OF BARE SOIL- NO RUNOFF 

17. TITLE VAN BAVEL - LASCANO 21 MAY 1979 

18. TITLE PROGRAM USED FOR USER'S GUIDE PREPARATION 

19. TITLE CONSERVB VARIATXON46 

20. 4*4* 

21. STORAGE T COM{ 25 ) , ITHET A ( 25 ) . DEPTH! 25 ) . CONDI 25 ) . HPOT ( 25) 

22. STORAGE A VCOND ( 2 5 ) . FLUX! 25 ) . PPO T ( 25 ) . 0 I ST ( 25 ) 

23. STORAGE AVKONO! 25 ) . FLOW ! 25 ) . KONO ! 25 ) . I TEMP! 25) , VHCAP I 25) , I VHCAO! 25) 


2A. / OlMENSiq'N VOL W! 25) . I VOLW ! 25 ) . NFLUX ( 25 ) , THET A ! 25 ) 

25. • / DIMENSION TEMP ! 25 ) . I VOLH! 25 ) . VOLH ! 25 ) . NFLOW! 25 ) 

26. / DIMENSION WINPUTU1.37) 

2?. /■ DIMENSION INF!38).RN!38).EVP! 38) .DRAINI38) 

28. / EQUIVALENCE ! VOLWl , VOLW! 1 ) ) , I I VOLW 1 , I VOLW ! 1 ) ) , I NFLU X I , NFLUX I 1 ) ) 

29. / EQUIVALENCE ! VOLHl , VOLH ! 1 ) ) . ! I VOLHl , I VOLH ! 1 ) ) , ! NFLOW 1 . NFLOW ! I )) 

30. FIXED I . J ,K , N) , Nl I , DNIIM, DNUM I . PW'.'MZ 

31. 4444444444444 * 44444444444444444 * 444444 * 4 * 44444444 ** 444*444444 

32. 444 * 

.33. 4*44 INITIAL SECTION 

34. 444* 

35. 4444 

36. INITIAL 

37. *44****4*4**44*4*44******444*44**44*»**4444*4**4**4******4444 

38. NOSORT 


39. 

4*44 1 ) 

OE 

FINITION OF PARAMETERS 

40. 

PARAMETER 


NL=13 

41 . 



NLL= NL *■ 1 

42. 

PARAMETER 


KONDS =4.2 

43. 

PARAMETER 


KONOW = 0.57 

44 . 

PARAMETER 


KCNOA = 0 .025 

45. 

PARAMETER 


VHCAPS= 1.925E06 

46. 

PARAMETER 


VHCAPW= 4.186E06 

47. 

PARAME TER 


SIGMA = 5.67E-08 

48 . 

PARAMETER 


20 = 0.01 

49. 

PARAMETER 


SATCON= 0.50E-06 

50 . 

PARAME TER 


PORSTY= 0.42 

51 . 

4444 2) 

DAILY COUNTERS 

52. 



DNUM = 1 

53. 



DNUM1= 2 

54. 



DNUM2= 3 

55. 

4*4* 3) 

READ INPUT DATA 

56 . 



DO 10 K= 1,37 

57. ' 



READ! 5.1000) IWINPUT! J.K) . J=1 . 1 1 ) 

58. 

1000 

FORMAT !F 5.0. 1X.F5.2,9!1X,F4.1)) 

59. 



IF! W I NPUT ! 1 , K) .EQ.O .0 ) CO TO 11 

60. 


10 

CONTINUE 
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61 . 
62. 

63. 

64. 
65 . 
66 . 

67. 

68 . 

69. 

70. 

71 . 

72 . 

73. 

74. 

75. 
76 . 
77. 

' 78. 

79. 

80. 
81 . 
82. 
83. 
84 . 
85. 
86 . 
•87. 
88. 
69. 
90. 
91 . 

92. 

93. 
94 . 
95. 
96 . 

97. 

98. 

99. 
100 . 

1 01 . 
102 . 
103. 

1 04. 
105. 

1 06 . 

1 07 . 
108. 

109. 

110 . 
111. 
112 . 

113. 

114. 

1 15. 
116., 
117. 

1 1 8 . 

1 19. 
120 . 
121 . 


11 CONTINUE 

4) SPECIFICATION OF THE GEOMETPY OF THE SYSTEM, 

INITIAL WATER CONTENT AND TEMPERATURE 

TABLE TCOMI 1-13) =0.01. 0.02. 0.02. 0.03. 0.03. 0.04, 0.05. 0.05.... 

0.10. 0.1 0.0. 15, 0.30. 0.30 
TABLE ITHETAI 1- 13 )=1340 .42 

TABLE t TEMPI 1-1 3) =3.57.4.92,6 .52. 8. 1 7, 9.66. 10.84, 1 1 .69, .. . 

12.08, 12. 18, 12.26, 12.69, 13,81 . 15, 22 

5) CALCULATION OF DISTANCE AND DEPTH 

OEPTHIl) = 0.54TCOMI1) 

DISTI 1) = DEPTH! 1) 

OO 20 1=2,NL 

OISTII) = 0,5*(TCOM( I-l )4TC0M( 1 ) ) 

OEPTH(I) =DEPTH(I-1) + 0 . 5* I TCOM ( I - 1 » +TC OMI I ) ) 

20 CONTINUE 

*♦** 6) CALCULATION OF INITIAL WATER AND TEMPERATURE CONDITIONS 

*♦** FLUX REFERS TO WATER 

FLOW REFERS TO HEAT 
IWATER = 0.0 
DO 30 I =1 ,NL 

IVHCAPin = VHCAPW41THETAI I )«■{ 1 .0-PORSTY)*VHCAPS 
IWATER =' IWATER TCOM ( 1 ) * ITHET A( 1 1 

IVOLH(I) = ITEMPI I )*TCOMt 1 )*IVHCAP( I ) i 

NFLUXI I ) = 0.0 
NFLOWI I > = 0.0 

IVOLWIIl = ITHETAI I )*TC0M( I ) 

30 CONTINUE 

♦ 444 

4444 HYDRAULIC CHARACTERISTICS OF FIRST HORIZON 0.0-0.20 M 

4444 volumetric WATER CONTENT V3. PREOCURE POTENT ! *.L 

4444 IN H WATER COLUMN 

FUNCTION TVSPl = < 0.030. -29000.00). ... 

I 0.050, -16000.00). ... 

( 0.070. -10000.00). ... 

I 0.090, -6000.00), , . 

I 0.110. -3000.00). ... 

I 0.130. — 1500.00). ... 

I 0.150. -700.00), ... 

I 0.170. -370.00), ... 

I 0.190. -160.00). ... 

I 0.210. -74.00) .... 

( 0.230. -45.00), ... 

I 0.250. -30.00). ... 

I 0.270. — 18.00). ... 

( 0.290. -11.00). ... 

( 0.310. -7.20). ... 

I 0.330. -3.60), ... 

I 0.350, — 1.20). ... 

I 0.370. — 0.60). ... 

( 0.390, -0.22). ... 

10.410. — 0.07).... 

( 0.420, 0,00), ... 

< 1.000. 0.00) 

4444 VOLUMETRIC WATER CONTENT VS. HYDRAULIC CONDUCTIVITY IN M/S 

FUNCTION TVSCl =' I 0.040. 0 . I 980 432E- 1 8 ) , ... 

( 0.060. 0.1867095E-1 7) . ... 

( 0.060. 0 .92 1 51 1 5E-1 7 ) . ... 

I 0.100. 0.3565G05E-16) . ... 

( 0.120. 0 . 14 1 1 493E- 15 ) , ... 

I 0.140. 0.6032844E-15) , ... 
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122 . 
123. 
124 . 

125. 

126. 

127. 

128. 

129. 

130. 
131 . 

132. 

133. 

134. 

135. 

136. 

137. 

138. 

139. 

140. 

141. 

142. 

143. 

144. 

145. 

146. 

147. 

148. 

149. 

150. 
151 . 

4 • 

153. 

1 54 . 
155.. 

156. 

157. 

158. 

159. 

160. 
161 . 
162. 

1 63 . 

164. 

165. 

1 66 . 
167. 

1 68 • 

1 69. 

170. 

171. 

172. 

1 73. 

1 74. 

175. 

176. 

177. 

178. 

179. 

1 60 . 
181. 

1 82 . 


( 0.160. 0 .2846457E- 14 1 , ... 

( 0.180. 0. 126331 7E-13) . ... 

( 0.200. 0.6312881E-13) . ... 

( 0.220. 0.3265092E-12) . ... 

( 0.240. 0. 1312614E-1 1 I . ... 

( 0.260. 0 .41 95633E- 1 1 ) . ... 

( 0.280. 0. 1 23991 4E-1 0) . ... 

( 0.300. 0.3590753E-10) . ... 

( 0.320. 0.9914647E-10I . ... 

( 0.340. 0.3065372E-09) . ... 

( 0.360. 0.1694214E-08) . ... 

( 0.380. 0.8728055E-08) . ... 

< 0.400. 0.5553067E-07 ) . ... 

( 0.420. 0.5000000E-06 ) . ... 

( 1.000. 0 .SOOOOOOE-06 ) 

HYDRAULIC CHARACTERISTICS OF SECOND HOR I ZON 0.20-1.20 M 
*♦** VOLUMETRIC WATER CONTENT VS. PRESSURE POTENTIAL 


FUNCTION TVSP2 


I 


0.020. 

-760.00) , 

• » • 

0.040. 

-340. 00) , 

• • • 

0.060, 

-130.00), 

» • • 

0.080. 

-60.00) . 

• • • 

0.100. 

-24.00 ) . 

• • • 

0. 120. 

-15.00). 

• • • 

0.140. 

-10.00) . 

» • » 

0.1 60. 

-7.40) , 

• • • 

0. 180. 

-5.40) , 

• # » 

0.200, 

-4.40). 

• » • 

0.220. 

-3.40) , 

• • • 

0.240. 

-2.80), 

• • • 

0. 260. 

-2.30) , 

• • • 

A AAA 

V 

-1.90. 

• • • 

0.300, 

-1.50), 

• • • 

0.320. 

-1.20) , 

• • • 

0.340. 

-0.94 ) , 

• • • 

0.360. 

-0.70), 

• • • 

0.380. 

-0.45) , 

V • • 

0.400. 

-0.27) . 

• • • 

0.420, 

-0.15). 

• » • 

0.440. 

-0.06 ) , 

• • • 

0.450, 

0.00). 

• • • 

1 .000. 

0.00 ) 



**** 


VOLUMETRIC WATER CONTENT VS. 


HYDRAULIC CONDUCTIVITY 


FUNCTION TVSC2 


0.030. 0 . 1200167E- 14 1 , 

0.050, 0. 1 599529E-1 3 ) , 

0.070. 0. I 51 6890E- 1 2 ) , 

0.090. 0 . I 06200 1 E- I 1 ) . 

O.llO, 0. 7476516E- I I 1 . 
0.130, 0 .34 70373E- 1 0 ) , 

0. 150. 0. 1 199559E-09) , 

0.170. 0.337535IE-09) . 

0.190. 0 .833201 7E-09 1 , 
0.210, 0. 1842094E-08) . 

0.230, 0 .376061 5E-08 ) , 

0.250, 0. 721 881 6E-08) , 
0.270, 0.1315320E-071 , 

0.290, 0.2300420E-07 ) , 

0.310, 0.3924749E-07 1 , 

0.330. 0.6606643E-07 ) , 

0.350. 0. 1 I 05985E-06 1 , 

0.370. O.I 068858C-06 ) , 
0.390, 0. 333631 9E-06 ) . 
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1 83 . 
184. 

1 85. 
186. 

187. 

188. 
189. 
1 90. 
191 . 

192. 

193. 

1 94 . 

195. 

196. 

197. 

198. 

199. 

200 . 
201 . 
202 . 

20 3. 

204. 

205. 

2 06. 

207. 

208. 

209. 

210 . 

21 1 . 
212 . 

213. 

214. 

215. 

216. 

217. 

218. 
219. 
220 . 
221 . 
222 . 

223. 

224. 

225. 

226. 
227 . 
22 8 . 
229. 

230 . 

231 . 

232 . 
233. 
234 . 

235. 

236. 

237. 

238. ^ 

239. ' 

240. 
24 1 . 

242. 

243. 


( 0.410. 0.671 1355E-06 1 . 
( 0.430, 0.1617800E-05) . 
( 0.450. 0.6000000E-05) . 
C 1.000. 0.6000000E-05) 


»*** 

SOIL TEMPERATURE VS. 

HEAT CONDUCTIVITY BY VAPOR IN W/IM4C) 

FUNCTI ON 

TEVSKO = 


-1 .000. 

0.02000) . 

• • • 




0.000. 

0.02470) . 

• • • 




10.000. 

0.04190 ) , 

• • • 




20.000. 

0. 07990) . 

• • • 




30.000. 

0.12600) . 

• • » 




40.000. 

0.24700) . 

• • • 




50.000. 

0.38100) , 

• • • 




60.000, 

0.65000) . 

• • • 




70.000, 

1. 17000) 


**** 

VOLUMETRIC 

WATER CONTENT OF FIRST LAYER VS. ALBEDO 

FUNCTION 

TIVSAL = 

! 

0.00. 0 . 

22 } t • • • 




( 

0.10. 0. 

22 ) • • • • 




! 

0.25. 0. 

1 7 ) • • • • 




! 

1.00. 0. 

17) 



**** 

**** 

**** 


7) TABLE OF THE GEOMETRY OF THE SYSTEM AND ITS INITIAL STATE 
NRITE(6. 1 100 ) 

1100 FORMAT! ‘0 I TCOM DEPTH ITHETA I TEMP*) 

00 40 1=1. NLL 

40 WRITE(6, 1200) I , TCOM U ) , DEPTH ( I ) , ITHET A( I ) , I T EMP ( I ) 
1200 FORMATdH .I2.4F10.5) 


******** 

**** 


******** *t ********* **tt********* ******* ****** ******** 


**** DYNAMIC SECTION 

**** 

V*** 

DYNAMIC 

************************************************************* 

NOSORT 

**** 

**** 1) DEFINITION OF TIME RELATED VARIABLES 

HTIME=TIME/3600. 

STIME=AMOO(HTIME.24. ) 

V=IMPULSI86400.,86400. ) 

IF(Y.LT.O.S) GO TO 22 • 

ONUM = ONUM + 1 

ONUMl = ONUMI ♦ 1 
ONUM2 = DNUM2 *■ 2 
22 CONTINUE 

XDNUM=FLOAT ( ONUM) 

JONUM = W INPUT (l.DNUM) 

**** 2) CALCULATION OF HYDRAULIC CHARACTERISTICS OF FIRST HORIZON 

DO 50 1=1.8 

THETA! I) = VOLW! I )/TCOM ! I ) 

CONO! I ) =AFGEN! T VSC 1 . THE T A ! I ) ) 

PPOT ! I I =AFGEN! TVSP 1 . THET A! I ) ) 

HPOT!I )=PPOT!I)-OEPTH(l) 

50 CONTINUE 

3) CALCULATION OF HYDRAULIC CHARACTERISTICS OF SECOND HORIZON 
DO 60 I=9.NL 

THETA! I) = VOLW! I >/TCOM I I ) 

CONO! I )=AFGEN! TVSC 2 . THE TA ! I ) ) 

PPOT ! I ) =AFGEN! T VSP2 . THET A ! I ) ) 

HPOT! 1 ) =PPOT { 1 ) -DEPTH! I ) 
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244 . 
245. 
246 . 
24 7. 

248. 

249. 

250. 
251 . 

252. 

253. 

254. 

255. 

256. 

257. 

258. 

259. 
260 . 
261 . 
262. 

263. 

264. 

265. 

266. 

267. 

268. 
269 . 
270. 
27 1 . 

272. 

273. 
274 . 

275. 

276. 

277. 

278. 
279 . 
280. 
281 . 
282. 
283. 
284 . 

285. 

286. 
287. 
286 . 
2 89 . 
29 0 . 

29 ) . 
292. 
293 . 
2 94 . 

295. 

296. 

297. 

298. 

299. 
30oL 

30 1 . 
30 2. 
303. 
30 4 . 


444* 4} 


4444 S) 


*444 6) 

4444 


44*4 7) 


4444 8) 

4444 


4444 9) 

4444 

4444 

4444 10) 
4444 11) 


60 CONTINUE 

CALCULATION OF THERMAL PROPERTIES 
DO 70 I=1«NL 

VHCAP(I) = VHCAPW4THET A( I ) 4 (1.0 - PORST Y ) 4VHCAPS 

TEMPI I) = VOLH< I ) /( VHCAPI I ) 4TC0M( I ) ) 

KdNDV = AFGEN ( TEVSK 0, TEMP ( I ) ) 

KONO(I) = (<I. -PORSTY) 4K0NDS40.4 4 THETA ( I ) 4K0NDW 4 ... 

(PORSTY - THETA I I ) ) 41 ,44( KONDA 4 KONDV))/ ... 

{ < 1 .-PORSTY )4 ,44THETA( I )4 ( PORST Y -T HET A ( I ) ) 41 .4 > 

70 CONTINUE 

CALCULATION OF AVERAGE CONDUCTIVITIES 

00 80 I = 2.NL 

AVKONO(I) = (TCOM( I-l ) 4TC0M( I ) )/( TCOM( I- 1 l/KONOt I -1 ) ... 

4 TCOMI I l/KONOI I ) } 

AVCONO( I.)=(CONO( I-l ) 4 TCOMI 1-1 )4CONDI I )4TC0M{ I ) )/. .. 
(TCOMI I-l )4TC0M( I ) ) 

80 CONTINUE 

SPECIFICATION OF BOUNDARY CONDITIONS FOR FLOW OF HEAT 
AND FLUX OF WATER 

FLUX(NLL)= CONDI NL) 

FLOWINLL) =(TEMP(NL)-ITEMP(NL) ) 4K0N0I NL ) /( TCOM I NL ) /2. ) 
CALCULAT ION. OF FLOW OF HEAT AND FLUX OF WATER 
DO 90 I = 2.NL 

FLOW! I ) = (TEMPII-1) - TEMP ( I ) ) 4 aVKOND ( I) /D I ST I I ) 

FLUX! I )=(HP0TI I- I l-HPOTI I ) )4AVCONOt I)/DISTI I ) 

90 CONTINUE 

USE OAILY RAINFALL DATA, AND CALCULATE DAILY 
RAINFALL DISTRIBUTION 

BEGIN = WINPUTI 9.0NOM) 

END = WINPUT(IO.ONUM) 

RT.T = WINPUT ( 1 I iONUM) 

RAIN=0.0 

IFIRFT.EQ.O.O) GO TO 33 
UPSLOP=I4.04RFT)/((END-8EGIN)442) 

OWSLOP=-UPSLOP 
M0PNT=(BEGIN4EN0 )/2.0 
H£IGHT=( 2.04RFr)/(ENO-BEGIN ( 

IFI ST IME .GE .beg IN. AND. ST I ME .LE.MDPNT) RAIN=. . . 

I UPSL0P4 { ST 1 ME-BEGIN > 1/3600000.0 
IFI STI ME .GT.MDPNT . AND. ST 1 ME . LE . END ) P A I N= . . . 
IOWSLOP4ISTIME-ENO) 1/3600000.0 
33 CONTINUE 

USE DAY LENGTH AND CALCULATE DAILY DISTRIBUTION OF 
GLOBAL RADI AT ION 

DL=WINPUT( 2.0NUM) 

OGR/ 86400 .41 .E06424 . /DL 4P I / 2 . = 4 36 . 33*DGR/DL 

GR=4 36. 334WI NPUTI 3. DNUM) /DL4S IN ( I ST I ME-1 2 . 4DL/2. ) . . . 
43.141/DL) 

IF (GR.LE.0.0) GR = 0.0 
CALCULATION OF ALBEDO 
T1 = THETAll) 

ALB =AFGENITIVSAL.T1 ) 

CALCULATION OF WINOSPEED AND BOUNDARY LAYER RESISTANCE 
IFI HT IME.LE . 12. ) SA=WI NPUTI 8, DNUM) 

IFIHT IME.LE.12. ) GO TO 44 

IF 1ST IME.LE. 12. )SA=WI NPUT I 8 . DNUM- 1 ) 4 I STI ME41 2 . 1/24 .4 . . . 

1 WINPUTI 8. ONUM)- WINPUTI 8. DNUM- 1 ) ) 

IFIST IME.LE.12. ) GO TO 44 

SA=WI NPUTI 8. DNUM) 4 ( ST I ME - 1 2 . ) /2 4 . 4 I W 1 NPUTI8 .DNUM4 1 ) - . . . 
WINPUTI 8, DNUM) ) 

44 CONTINUE 
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305. 

3C6. 

307. 

3Ce. 

309. 

310. 

311. 

312. 

313. 
314 . 
31 S. 
316 . 

31 7. 
316. 

319. 

320. 
321 . 

322. 

323. 

324. 

325. 

326. 

327. 

32 8. 

329. 

330. 
331 . 
332. 
^ ^ 

334. 

335. 

336. 

337. 
336. 

339. 

340. 

341. 
342 . 
343. 
344 . 

345. 

346. 
347 . 
346. 
3 49 . 
350. 
351 . 
352. 

353 . 

354 . 

355. 

356. 

357. 

358. 

359. ' 

360. 
361 . 

362. 

363. 

364. 

365. 


**** 


**** 


* 44 * 

4444 


4444 


4444 


4444 


RA e (ALOG(2.0/ZO|442.0)/(0.164SA) 

12) CALCULATION OF DEWPOINT TEMPERATURE AND ABSOLUTE HUMIDITY 

0PMAX=WINPUT(6.DNUM) 

DPMIN=WI NPUt (7 , ONUM) 

DPTC=DPMIN4(DPMAX-0PMIN)4 (STIME-5. )/l 0. 

1F(ST IME.GT . 15. ) OPM I N=W I NPUTt 7, 0NUM4I ) 

1F< STIME.GT.15. )DPTC=DPMAX-(0PMAX-0PMIN»4 (STIME-15. )/l 4. 
IF(ST IME.LT .5. AND .0NUM.GE.2. 1 DPH A X=W I NPUT ( 6, ONUM- 1 ) 

XF( STIME.LT.S. )0PTC=DPMAX-(DPMAX-0PMIN)4(STIME49 . )/l4 . 

HA s 1.3234EXPI 17.274DPTC/(237.34DPTC))/(273.I64DPTC) 

13) CALCULATION OF TEMPERATURE OF THE AIR AND SH OF THE AIR 

TAMAX=W INPUT (4 .ONUM) 

TAMIN=WINPUT (5. ONUM) 

TAC=TAMIN4(TAMAX-TAMIN) 4{ STI ME-5. )/10 . 

IF ( ST IME.GT . 15. )T AMI N=W INPUT! 5. DNUM+ 1 ) 

IFfSTIME.GT.15. ) T AC = TAMAX- ( T AM AX-T AM I N ) 4 ( ST 1 ME-1 5 . )/l 4 . 
IF(ST IME.LT .5.AN0 .ONUM.GE .2. ) TAMA X=WI NPUT ( 4 .ONUM- 1 ) 

IFf ST IME.LT. 5. )TAC = TAMAX-( TAM AX-T AMIN)4(STIME49.)/14. 

TAK=TAC4273.16 

SH=( 1 154 .84303. 16)/ ( TAX ) 

14) CALCULATION OF SKY IRRADIANCE 

SKL = ( SI GMA4TAK444 )4( 0.60540 .0394SQRT< 1410. *HA)) 

15) IMPLICIT CALCULATION OF THE SOIL SURFACE TEMPERATURE 

TSC = IMPL (TAC.0.01 ,FTSC ) 

A = (TSC - TA04SH/RA 

HO a 1.3234EXP! 17.274TSC /(237.34TSC ) ) / ( 273 . I 64TSC ) 
HO=H04EXP ( PPOT ( 1 ) /( 46 . 974 ( T SC 42 73 . 1 6 ) ) ) 

EV a(H0 - HA)/(RA41000. ) 

S a GR4(1. - ALB) 4 SKL - S1GMA4ITSC 4 273.16)444 ... 

- A -LH4EV 

FTSC a TEMP! I) 4 S40EPTH ! 1 ) /KOND! 1 ) 

16) CALCULATION OF EVAPORATION AND NET RADIATION 

FLOWIl) a (TSC - TEMP! 1 ) ) 4K3N0! 1 ) /D1 ST( 1 ) 

HS a 1 .3234EXP! 17.274TSC /!237.34TSC ) ) /( 273 . 1 64T SC ) 

HSO = HS 

HS=HS*EXP!PPOT! I )/!46.974(TSC4273.l6) ) ) 

RHS a HS/HSO 

EVAP a (HS - HA)/!RA41000. ) 

LH=2. 49463E09-2.247E064TSC 

NR a FLOW!l) 4 !TSC - TA04SH/RA 4 LH4EVAP 

17) CALCULATION OF DETAIN, INFILTRATION, AND INCAP 

DETAIN a INTGRL !0.0. RAIN-INFILT) 

INCAP a (O.-HPOT! 1 ) )*0.54(SATCON4COND! 1 ) ) / OISTIl) 

IF (RAIN. GT. 0.0) GO TO 55 
IF (DETAIN, LE. 0.0) GO TO 66 
INFILTaINCAP 
GO TO 77 
66 CONTINUE 

DETAIN a 0.0 
INFILTaO.O 
GO TO 77 
55 CONTINUE 

INFILT a INCAP 

IF (RAIN. LT. INCAP. AND. DETAIN. LE . 0 . ) I NF I LTap A I N 
77 CONTINUE 

18) CALCULATION OF NET FLOW OF HEAT AND NET FLUX OF WATER 

FLUX( 1 )=I NFILT-EVAP 
DO 100 I a l.NL 

NFLOW(I) a FLOWIl) - FL0W(l4l) 

NFLUX! I )aFLUX( I) -FLUX ( I 41) 

100 CONTINUE 
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366. 

367. 
366. 

369. 

370. 
371 . 
37 2. 

373. 

374. 
375 . 

376. 

377. 

378. 

379. 

380. 
381 . 

382. 

383. 

384. 

385. 

386. 

387. 

388. 

389. 

390. 

391. 

392. 
39 3. 
394 . 
395. 
396 . 
397. 
39 8. 

399. 

400. 
401 . 
402. 

403 . 

404 , 

405. 

406. 

407. 

408. 

409. 
410 . 

411. 

412. 

413. 

414. 

415. 
4 16. 

417. 

418. 
4 19. 
420. 

421 . 

422 . 
423. 
424 . 

425. 

426. 


19) 

44>«4 


20 ) 

**** 


**** 21 ) 


I 


22 ) 


INTEGRATION OF VOLUMETRIC HEAT CONTENT AND VOLUMETRIC 
WATER CONTENT 

VOL HI =INTCRL( IVOLHl .NFLOWI ,13) 

V0LW1 = 1 NTGRLI IVOLWl . NFLUXl , I 3) 

CALCULATION OF CUMULATIVE RAIN. INFILTRATION, EVAPORATION 
AND DRAINAGE 


CUMRN = INTGRL 

! 0 

. 0 . 

RAIN) 

CUMINF = INTGRL 

! 

0.0 

. INFILT ) 

CUMEVP = INTGRL 

1 

0.0 

. EVAP ) 

CUMORN = INTGRL 

! 

0.0 

, FLUX! NLL ) ) 


CALCULATION OF TOTAL WATER FOR DIFFERENT LAYERS 
CUMWTR = 0.0 
CUMWTl = 0.0 
CUMWT2 = 0.0 
CUMWT3 = 0.0 
OO no I=1.NL 
CUMWTR= CUMWTR ♦ VOLW(I) 
no CONTINUE 

DO 120 I=6.NL 

CUMWTl = CUMWTl ♦ VOLW(I) 

120 CONTINUE 

OO 130 I=9.NL 

CUMWT2 = CUMWT2 ♦ VOLWt I ) 

130 CONTINUE 

OO 140 I=12,NL 

CUMWT3 = CUMWT3 VOLW( I ) 

140 CONTINUE 

CALCULATION OF DAILY TOTALS 

ZBHJS=IMPULS(864 00.,864 00. ) 

IFIZBHJS.lt. 0.5) GO TO 88 
; Nr { 0NUf4 J ) =CUM1 NF 
RN(DNUM-I )=CUMRN 
EVP(ONUM-l ) sCUMcVP 
DRAIN! ONUM- 1 )=CUMORN 
OINF=CUMINF-INF (ONUM-2) 

DRN=CUMRN-RN( DNUM-2 ) 

OEVP=CUMEVP-EVP( ONUM-2 1 
OORN=CUMDRN-OP Al N( ONUM-2) 




23) 

24) 


S 


$ 


88' CONTINUE 

CHECK WATER BALANCE 

BALANS = CUMWTR - IWATER - CUMINF CUMEVP «• CUMORN 
OUTPUT OF DESIRED VARIABLES 
Z=IMPULS(2 8800.0,864 00. ) 

IF(Z.LT.O.S) GO TO 99 
222 CONTINUE 

WR1TE(6 , 1 300) WINPUT! 1, DNUM) .TIME, XONUM.STIME 
1300 FORMAT!* JULIAN DAY NUMBER = '.FA.O.* TIME = '.FlO.l 
• XDNUM = *,F11.0,* STIME = '.FT.A) 

WRl TE!6, 1400) 

1400 FORMAT!'0 I • . 5X . • DEPTH • , 1 0 X , • THE T A • , 1 1 X , • PPOT • , 1 1 X , 
•FLUX* ,9X . *NET FLUX* , 1 OX, • TEMP* ) 

DO 150 1=1, NL 

150 WRI TE!6 , 1500 ) I . DEPTH ! I ) , THET A ! I ) , PPG T ! I ) , FLU X ! I ) , 
NFLUX! I ) , TEMP! I ) 

1500 FORMAT! I3.6E15.4) 

WRl TE!6 . 1600) 

1600 FORMAT! IH ////) 

GO TO 111 
99 CONTINUE 

Z2=I MPULS !576 00. ,86400.) 

IF! Z2.LT. 0.5) GO TO 111 
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427. 

426. 

429. 

430. 
431 . 
432. 

433 . 

434 . 

435. 

436. 

437. 

438. 

439. 
440 . 
44 1 . 
442. 
442. 
444 . 
445. 
446 . 

447. 

448. 

449. 

450 . 

451 . 

452. 

453. 

454. 

455. 

456. 
457 . 
458. 
’459. 
460. 

461 . 

462 . 
463. 
464 . 

465. 

466. 


GO TO 222 
111 CONTINUE 


**** 

***♦ 

4444 


E TERMINAL SECTION 


terminal 


TIMER F NTIM=259200.0,PROEL=86400.0.0ELT=100.0 
PRINT JONUM. XONUM, ORN. OINH. OEVP. ODRN. CUMWTR , . . . 

CUMWTl . CUMMT2. CUMWT3. BALANS. CUHRN. CUMDRN . CUM INF . CUMEVP , . . . 
FLOW! 14) 

METHOD TRAP2 
END 
4 
4 
4 


* WEATHER INPUT DATA, STORED IN ARRAY W1 NPUT ( 1 1 , 37 » 


4 JNM 
INPUT 

OL 

OGR 

TMAX 

TMJN 

DMAX DM1 N 

SA 

BEGIN 

END 

RFT 

1 . 

10.70 

10.0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

O 

• 

o 

0.0 

2. 

.10. 70 

1 0. 0 

16.0 

0.5 

3.0 -2.5 

3. 0 

9.0 

12.0 

25.0 

3. 

10.70 

10.0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0.0 

4. 

10.70 

10. 0 

16. 0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0.0 

5 . 

10.70 

10.0 

16.0 

‘ 0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0.0 

6. 

10.70 

10.0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0.0 

7. 

10.70 

10.0 

16.0 

o‘.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0 .0 

8. 

10. 70 

1 0.0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0 .0 

9. 

10.70 

10.0 

1 6.0 

0.5 

3.0 -2.5 

3. 0 

0.0 

0.0 

0 .0 

10. 

10.70 

10.0 

16.0 

0 .5 

3.0 -2.5 

3. 0 

0.0 

0.0 

0.0 

1 1 . 

10.70 

1 0. 0 

1 6. 0 

0.5 

3.0 -2.5 

3.0 

0.0 

o 

• 

o 

0.0 

12, 

10.70 

10.0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0 • 0 

13. 

10. 70 

1 0. 0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0.0 

14. 

10.70 

10.0 

16.0 

0.5 

3.0' -2.5 

3.0 

0.0 

0. 0 

0.0 

15. 

10.70 

10.0 

1 6.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0 .0 

0. 10.70 
ENOINPUT 

STOP 

END JOB 
/4END 

10.0 

16.0 

0.5 

3.0 -2.5 

3.0 

0.0 

0.0 

0.0 
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APPENDIX C - CALCULATION OF UNSATURATHD HYDRAULIC CONDUCTIVITY 
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CALCULATION OF UNSATURATED HYDRAULIC CONDUCTIVITY 


PURPOSE : Uiis WATFIV program is designed to calculate the unsaturated 

hydraulic conductivity as a function of water content using Jackson's 
method (Jackson, 1972: Bouwer and Jackson, 1974). It also states the 
relationship between water content vs. pressure potential and water 
content vs. hydraulic conductivity in the form required for a CSMP III 
FUNCTION statement. Note, that the synfcols used for variables in this 
program are different from those used in the general model. 

INTRODUCTION : In order to use this program, three parameters must be 

known. These are: 

1. the saturated hydraulic conductivity of the soil, (CCWDCl)), in 

ni/s, 

2. The relationship between pressure potential, (PP(I)) , in m water 
column and volumetric water content, (TH(I)), and 

3. the maximum value of water content, (1H(1)). 

The first two parameters are obtained experimentally and the third one 
is estimated by calculating the porosity of the soil, i.e. 

TH(1) = porosity = 1 - DB/PD 

where DB is the dry bulk density of the soil in g/cm^, and PD is the 
particle density of the soil, usually taken as 2.65 g/cm^. 

It should be noted that Jackson's method requires detailed 
information on the relationship between pressure potential vs . water 
content for high values of pressure potential, i.e.> - 0.50 m. 

Experience has shown that this information has a critical effect on the 
results . 
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PROCEDURE : This section outlines the procedure that should be followed 

to obtain the values of pressure potential for different values of water 
content from experimental data. 

1. Plot the relationship between pressure potential vs. volumetric 
water content for values of pressure potential > - 5.0 m on a linear 
scale. To the value of 0.0 pressure potential assign the maximum 
water content (TH(1)). A semi- log graph of pressure potential vs. water 
content is used for the remaining values of pressure potential. 

2. Select an increment for the volumetric water content (DELTH) ,• 
i.e. 0.01, 0.02. 

3. Calculate the number of values of water content (M) for which the 
hydraulic conductivity is to be calculated. 

M ° 1H(MAXIMUNf) - TH(MINIMUM) i n 
DELTH 

4. From the graphs, read the values of pressure potential that 
correspond to the midpoint of ekch equal increment of volumetric 
water content. This is illustrated with an exanple. 


EXAMPLE: 

DELTH = 0.02, THCl) 

= 0.42 


WATER CONTENT 
INTERVAL 

MIDPOINT 

PP(I) c 

0.42 

- 0.40 

0.41 

0.07 

0.40 

- 0.38 

0.39 

0.22 

0.38 

- 0.36 

0.37 

0.60 

0.04 

- 0.02 

0.03 

29000.00 
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I NPUT TO THE WATFIV ITOGRAM: The following input cards must be specified, 

in the sequence given, after the //$DATA card. 


DATA CARD NO. 

IDENTIFICATION 

RESTRICTIONS 

1. 

TITLE 

up to 80 alphabetic 
or numeric charaters 

2. 

TH(1) 

Real variable 

3. 

DELTH 

Real variable 

4. 

M 

Integer 

5. 

PP(1) 

Positive real 
variables 

1+4 

PP(I) 


* 

M+4 

PP(M) 


M+5 

COND(l) 

Real variable 


OUTPUT ; The output of the program is in three parts. The first 
part lists the title and the input data. Note that in the first 
table the values of PP(I) , LOGIO PP(I) , and RH(I) correspond to 
THETAM(I), the midpoint water content. The second part lists the 
results. Note that in the second table the calculated value of hy- 
draulic conductivity, (CX)NDCI))» corresponds to TIETA(I) . The 
maximum hydraulic conductivity (CXDND(l)) should correspond to the 
maximum value of water content (TH(1)). The third part is a table of 
volumetric water content vs. pressure potential, and volumetric water 
content vs. hydraulic conductivity, in the form required for a CSMP 
FUNCTION statement. 
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1 . 
2 . 
). 

4. 

5. 
5. 

7. 

8. 

•). 

10 . 

11 . 

12 . 

13 . 

14. 

15. 
15. 

17. 

18. 

19. 

20 . 
21 . 
22 . 

23. 

24. 

25. 

26. 

27. 

28. 

2 9. 

30. 

31. 

32. 

33. 

34. 

35. 

36. 

37. 

38. 

39. 

40. 

4 1. 

42. 

43. 

4 4 . 
45. 

4 6. 
47. 

4 8. 

49. 

50. 

51. 

52. 

53. 

54. 

55. 

56. 
57.. 
59. 

59. 

60. 


//JACKSOB J08 (8042, 402E, ♦02,001, ?L) LASCAHO 

/♦WATFIV 

C 

0m 

w 

C CALCULATE HYDPAIILIC COUDHCTI VITY 

C JACKSON 1972 SSAP 36 ; 380-302 

t 


THB HSCESSAPY INPUT FOR THIS CALCULATION IS GIVEN BELOW 
THE EXPLANATION OF INPUT VARIABLES IS GIVEN IN THE PBOGRAH 


C 

C 

c 


C 

r» 

w 

c 

c 

c 

c 


c 

e 


DATA CAPO 

2 

3 

4 

5 


IDENTIFICATION 
. . ..TITLE 
...,TH(1) 

. . . .DSLTH 

. . . . n 

....pp(i) 


I *4. 


PP(I) 


H*4 


PP(Jl) 

COND(1) 


REAL NOai, N01.1, LCOHB, LPP 
INTEGER TITLE(OO) 

DIHENSION TH(701, PP(70j, ?ELCON(70), C04D(70), LCOND(70), ?H(70J 
DIMENSION LPP (70) ,TIIETA.'1 (70) ,THX(70) ,CONDX (70) ,PPX (70) , THY (70) 

INPUT DATA 


C READ TITLE (NAME OF SOIL, DEPTH ,OTII EP. INFOPMATION) 

READ{5, 1 1) TITLE 
11 FOPHAT («0 A1) 

C MAXinOn VALUE or water content 

READ.Tli (1) 

Z DEEINE STpf’ STZE (INCREMENT) 

READ, PELTH 

C M IS THE, TOTAL NUMBER O? VALUES OF WATER COMENT FOR WHICH THE 

Z HYDRAULIC CONDUCTIVITY IS CALCULATED 

: H CAN BE CALCULATED FROM THE FOLLOWING EXPRESSION: 

C M = (TH (MA7IMUM) -TH (HINI-UH) ) /DELTH ♦ 1.0 

READ, n 

C PPOGRAn 


DO 10 I = 1,H 

TH(I) s TH(1) - (I-1)*nElTH 
THETAM(T) = TH ( T) - (DELTH/2. 0) 

READ PRESSU.EE POTENTIAL VALUES IN M, CORRESPONDING TO THEIR 
RESPF.CTIVE VALUES FOR THE CENTER OF EACH INCREMENT THETA. 
READ,PP(T) 

PP(I) = -PP(I) 

RH(I) = EXP(PP(I) /14091.) 
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61. LPP(I) = ALOG10 (-PP(rj) 

52. 10 CONTINUP 

53. DEHOn = 0. 

50 . DO 20 J = 1 , H 

55. DENOMJ = (2*J-1) /(PP (.1) **2) 

55. DEMON = DEMON ♦ DEKCMJ 

57. 20 COMTINOE 

53. r CIVS NAXINiJN CONDUCTIVITY IN H/S (SAT. CONDUCTIVITY). 

59. PEAD, COMD(1) 

70. C 9ECOPD OF INPUT DATA 

71. WRITE(6,850) TITLE 

72. 850 FOFHAT {' 1 • ,//,20X, 'TITLE: *,80A1,//) 

73. MPITE (6, 900) 

74. 900 FOSRATC-' ,//, USX, 'INPUT DATA* ,//) 

75. MPITE(5, 1000) TH(1) 

75. 1C00 FOPNAT , 20X, ' HAXIMUH VALUE OF WATER CONTENT =»,F0.2) 

77. . WRITE(6, 1 100) DSLTH 

78. 1100 FORHATC-' ,20X, 'INCREHENT OF MATER CONTENT =',F5.3) 

79. WPITE(5, 1200) H 

30. 1200 FORNAT('-',20X, 'NUNBER OF VALUES OF HATER CONTENT, FOP WHICH X IS 

31. ICALCULATSD = *,I3) 

32. WRITE (6, 1 300) COND(1) 

33. 1 303 FOPNAT('-',20X,'NAXIM1IN CONDUCTIVITY (SAT. COND.) =',E14.5 

84. S,2X,'N/S') ; ■ > 

85. WRITE (6, 1400) 

36. 1400 FORMAT (• O' , 20X, ' THE VALUES OF THETAN APE THE MIDPOINTS',/, 

37. S21X,'OF EACH INCREMENT TO WHICH THE VALUES OF PP(I)',/, 


38. 


S21X, 'CORRESPOND ACCORDING TO EXPERINENTAL 

DATA') 


.39. 


WR.ITE(6, 1450) 



93. 

14 50 

FORNAr('0',20X, •TH*' lOGIO VALUE OF PP(I) 

AND RH (I) 

CORRESPOND' , 

9 1 . 


S/,21X,'TO THE VALUE OF THETAN ( I) ’ ,////) 



92. 


WPITE(6, 1500) 



93. 

1500 

FORHAT C-',20X, ' THFTA (I) ' ,2X, 'THETAN (I) ' , 

4X, ' PP (I) 

(K) ',2X, 

94. 


S'LOGIO P?(I)',2X, 'FH(I)') 



95. 


DO 50 1=1,. 1 



96. 


WRITE (6, 1500) Ttl (I) , THETAN (I) , °P (I) ,LPP(I) 

,PH(1) . 


97. 

1600 

.FORMAT ( ' ' ,23X, F4. 2,5X, F4. 2,5x,F9. 2, 6X,F6 

.2,2X,F6. 

2) 

93. 

50 

CONTINUE 



99 . 


DO 30 I = 1,H 



100. 


NOMI = 0. 




131. DO 40 J = I,N 


102. 


NOHJ = (2*.1 ♦ 1 - 2*1) /(PP (J) *♦ 2) 

133. 


NOr.I = NON’’ ♦ NOMJ 

194. 

40 

CONTINU'" 

135. 


HF.LCON(I) = (TH (I) /TH (1) ) * (NONI /DEMON) 

196. 


CO'.'DID = P.ELCON (I) =CCND ( 1) 

197. 


LCOND(I) = ALOG 10 (COND (I) ) 

193. 

33 

CONTINUE 

199. 


WRITE (6, 1700) 

119. 

1700 

FOP.IAr { • 1 ',//////, 4 5X, 'R'^'S'JLTS' ,//) 

111. 


WRITE (5 , 1 300) 

112. 

1800 

FOP NAT (' - ' , 20X, ' THETA (1) ' ,4 X, t reLCON (I) ' , 4X, ' COND (I) 

113. 


S'LOGIO CONO (I) ') 

114. 


DO 50 1=1 

115. 


MR:TE(6,1900) TH(I) , RELCON (I) , CON D (I) ,LCONO(I) 

115. 

1900 

FOP NAT (' ',21X,F6. 4,4X,S10.4, 3X,E10.4,4X,F6.2) 

117. 

60 

CONTINUE 

118. 


WRITE (6,2000) 

119. 

2000 

FORMAT(i:i1) 

120. 


K = 0 

121 . 


DO 90 J = 1 , N 
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122 . 

123. 

124. 

125. 

126. 
127. 
123. 
129. 

no. 

131. 

132. 

133. 

134 . 

135. 

136. 

137. 

138. 

139. 

140. 

141. 

142. 

143. 
14 4. 

145. 

146. 

147. 
143. 

149. 

150. 

151. 

152. 

153. 

154. 

155. 

156. 

157. 

158. 

159. 

160. 
161. 
162. 

163. 

164 . 

165. 

166. 
167. 
158. 

169. 

170 . 
171 . 


KK = :i-!1 

THX (J) =THF.TH?1<KK) 

CONOX (J) =COSO (KK) 

THY (J) =?H(KK) 

PPX (JJ =PP(KK) 

K = K*1 

9C CONTIHUB 

W»ITR(6, 2100) 

. 210C FOPHAT(llll) 

no 100 j^i.n 

3«»ITE (fi, 2200) THX (J) ,PPX (.1) 

2200 FOPMATI* ' , 20X, • ( * » • 5 » * » * » n 0 . 2 , ' ) , ...’) 

100 Cm.’TIMtIE 

Wr:TB(6,2300) 

2300 ?0n.1A?(1H1) 

DO 110 .7=i,n 

W»ITE(6,2400) THY (J) ,COKDX (J) 

2400 FOPIIATC •,20X, • (• «F5.3,‘ ,* »E14.7,») f 
110 CONTINUE 

W»ITE(6,2500) 

2500 FOB.5AT (1!I1) 

STOP 

END 

//$DATA 

NORWOOD SILT LO.M, AGRONOMY FARM, 0.00 - 0. 20 « DEPTH 
0.42 
0.020 
20 

0.07 
0 . 22 
0. 60 
1. 20 
3.60 
7.20 
11.0 
1B.0 

30.0 

45.0 

74.0 

160.0 

370.0 

700.0 

1500.0 

3000.0 

6000.0 

10000.0 

16000.0 

29000.0 
5.0E-07 
/♦END 
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ITLi: 


NORWOOD Slir LO/iM, AoKONOMY FARM, 0.00 


0.20 H DF.PTH 


ISrOT DATA 


.MAXIMUM VALUE OF WATER CONTENT =0.42 


INCREMENT OF WATER CONTENT =0.020 


NUMBER OF VALUES OF WATER CONTENT, FOR WHICH K IS CALCULATED = 


.MAXIMUM CONDUCTIV: 7 / (SAT. COND.) = 0.50000E-06 ^:/S 

THE VALUES OF THETAM ARE THE MIDPOINTS 
OF EACH INCREMENT TO WlilCH T!!S VALUES OF PP(I) 
COPRES?OND ACCORiilKG 70 EXPERIMENTAL DATA 

THE LOG10 VALUE 0.^ PP (I) AND ?H(I) CORRESPOND 
TO ThS VALUE OF THETAK(I) 


STAdJ 

THETAM (I) 

??(I) (M) 

LOG10 PP(I) 

PH (I) 

0.42 

0.4 1 

-0.07 

-1.15 

1. 00 

0.40 

0.39 

■ -0.22 

-0.66 

1 . 00 

0.38 

0. 37 

-0.60 

-0.22 

1. 00 

C. 3G 

0.35 

-1.20 

0.08 

1 . 00 

0. 34 

0.33 

-3. 60 

C. 56 

1 . 00 

0. 32 

0.31 

-7.2C 

0. 96 

1. 00 

0.30 

0.29 

-11.00 

1 . 04 

1. 00 

0.20 

0.27 

- 18.00 

1.26 

1. GO 

0.26 

0.25 

-30.00 

1.48 

1.00 

0. 24 

0.23 

-45.00 

1.65 

1. 00 

0,22 

0.21 

-74.00 

1,87 

0. 99 

O.TO 

0.19 

- 160.00 

2. 20 

0. 99 

0. 18 

0. 1 7 

-370.00 

2. 57 

0. 97 

0. 16 

0.15 

-700.00 

2.85 

0. 95 

0. 14 

0.13 

-1500.00 

3. 18 

0. 90 

0.12 

0- 1 1 

-3C0G.00 

3.4 8 

0.01 

0. 10 

0.09 

-6000.00 

3.78 

0. 6B 

0. 08 

0.07 

-10000. CO 

4.00 

0. U9 

0.06 

0.05 

-IbOOO.OO 

4. 20 

0. 32 

0.04 

O.OJ 

-29000.00 

4.46 

0. 13 
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RESULTS 


THETA (:) 

RELCOli <I) 

0.tl200 

0. lOOOh 01 

0.4000 

0. 1111L cJ 

0. 3000 

0. 1740E-01 

0. 3(>G0 

0. 3388E-02 

0. 3400 

0.61J1£:-(>J 

C. 3200 

0 . 1 '39 j I-;- j j 

0. 3000 

0 • 718^ L- tiU 

0.2000 

0. 248CL;-C* 

0.2600 

0.8391E-05 

0.2400 

0.2625E-03 

0.2200 

0.6530E-06 

0.2000 

0 . 126 j£- 0 o 

0. 1800 

0. 2527!>J7 

0. 1600 

0.569JL-C:3 

0. 1400 

0. 1207L-3'J 

0. 1200 

0. 282jE-03 

0. 1000 

0. 71J12-10 

C* OdOO 

C. ir.43T VO 

0.0600 

0.3734S-11 

0.0400 

0.3S61E-12 


COND(I) LOT 1C COSP(l) 


0.6000E-06 

-6. 39 

0.55:OE-07 

-7.26 

0. 8728E-03 

-G.06 

C. 1h Tur-Crt 

-3.77 

0. 30<^.^E-0'J 

-9.51 

U.431L.F.-10 

-10.00 

0.35012-10 

-10.44 

0. 124CE-10 

-1C. 91 

C.UISRE-II 

-11.38 

0. 13132-11 

-11.88 

0. 32G5E-12 

- 12. 40 

0.6J13E-13 

-13.20 

0. 12f 32-13 

-13.9) 

0.28462-14 

-14 .55 

0.6033E-15 

-15.22 

0. 1411E-1S 

-15. B5 

0. J566E-16 

-16.43 

n. 5215E-17 

-17.04 

0. 1H67E-17 

-17. 73 

0. 1900S-13 

-18.70 
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( 0.030, 
( 0,050, 
( 0.070, 
( 0.090, 

( 0. no, 

( 0.130, 
( 0.150, 
( 0. 170, 
{ 0.190, 
( 0 . 210 , 
( 0.230, 
( 0.250, 
( 0.270, 
( 0.290, 
( 0,310, 
( 0.330, 
( 0.360, 
( 0.370, 
( 0.390, 
( O.tilO, 


-29000. GO) , 
-16000. ut) , 
-10000.00) , 
-600C. OU) , 
-3000.00) , 
-1500.00) , 
-700. OJ) , 
-370.06) , 
-160.00) , 
-7U.0G) , 
- 45. 00) , 
-30.00) , 
-10.00) , 
-11 .00) , 
-7.20 . 
-3.60) , 
-1.20) , 
-0.60) , 
-0. 22) , 
-0.07) , 
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( o.ouo, 
{ 0.060, 
( 0.C80, 
( 0 . 100 , 
{ 0 . 120 , 
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